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ABSTRACT 

This paper describes a comprehensive measurement model for the error budget of pulse arrival times. 
Applications include forecasting timing errors for optimization of long-duration timing campaigns and 
identifying mitigation methods that can reduce the errors. The model takes into account many end- 
to-end effects from the neutron star to the telescope, with emphasis on plasma propagation effects 
(particularly interstellar scattering) , which are stochastic in time and have diverse dependences on radio 
frequency. To reduce their contribution, timing measurements can be made over a range of frequencies 
that depends on a variety of pulsar and instrumentation-dependent factors that we identify. A salient 
trend for high signal-to-noise measurements of millisecond pulsars is that time-of-arrival precision is 
limited cither by irreducible interstellar scattering or by pulse-phase jitter caused by variable emission 
within pulsar magnetospheres. Forecasts that ignore these contributions are likely to be overoptimistic. 
A cap on timing errors implies that pulsars must be confined to low dispersion measures and observed 
at high frequencies. Use of wider bandwidths that increase signal-to-noise ratios will degrade timing 
precision, after a point, if nondispersive chromatic effects are not mitigated. The allowable region in 
the dispersion measure-frequency plane depends on how chromatic timing perturbations are addressed 
in wideband measurements. Without mitigation, observations at 1.4 GHz or 5 GHz are restricted to 
DM < 30 and < 100 pc cm -3 , respectively. With aggressive mitigation of interstellar scattering and use 
of large telescopes to provide adequate sensitivity at high frequencies (e.g. Arecibo, FAST, phase 1 of the 
SKA, and the SKA), pulsars with DMs up to 500 pc cm -3 can be used in precision timing applications. 
We analyze several methods that involve wide-band fitting to arrival times at a given epoch prior to multi- 
epoch fitting. While the terms of greatest current astrophysical interest are achromatic (e.g. orbital and 
gravitational wave perturbations), their measurement may ultimately be limited by similarly achromatic 
stochasticity in a pulsar's spin rate. 



1. INTRODUCTION 



Timing of pulsars commenced with their original discovery over 40 years ago and has reaped enormous 
benefits in enabling fundam ental tests of general rela tivity ()Kramer et al.|[2006[ ). demonstrating the exis- 
tence of gravitational w aves (|Tavlor fc Weisberg1ll982ft , determining masses and other funda mental aspects 
of neutron stars fNSsI ([Freire et~l2009ft . and modeling of the interstellar medium (ISM) (|Cordes et al.l 
119901: iRamachandran et all I2OO60 . The precision of time-of-arrival (TOA) measurements has improved 
progressively, allowing new areas of study to be targetted, such as higher-order relativist ic corrections to 
orbits and the timing perturbation from low-fr equency (nano-Hertz) gravitational waves (lDetweilerlll97 
Hclli ngs fc Downslll983HRomani fc Tavlorl 19831 ). The discovery of millisecond pulsars fMSPsi (|Backer et al 



1982ft opened the door to placing meaningful limits on or making detec tions of cosmological g ravitational- 
wave (GW) background s (|Bertotti et al.lll983tlCordes fc Stinebrindll984l: [Foster fc BackeJl99CtlStinebring etaH 



1990; iKaspi et al.lll994ft . Projections based on improved timing met hods and discovery of new MSPs now 
makes the detection of nano-Hz GWs p lausible in the next decade (ICordes et al.l 12004 iJenet et al. 2004; 
Kramer et alj|2004 Uenet et alj |2005allbt iLee et alj|200§ Uenet et al.ll2009HHobbs et alJl2009HVerbiest etldl 



2009| iHobbs et al 



2010). 



This paper is motivated by the general goal of maximizing timing precision but especially for GW detection. 
The sensitivity of a pulsar-timing array (PTA) increases nonlinearl y with data span length but likely requires 
5 to 10 years to make detections of plausible GW backgrounds l|Anholm et al.ll2~009t Ivan Haasteren et al.l 
2009). The relevant GW frequencies are therefore in the 10 -8 — 10 -9 Hz range. Candidate signals in- 
clude stochastic backgr ounds from coalescences of supermassive black-hole binaries (jjaffc & Backer 2003; 
ISesana fc Vecchid 20101) in galaxy mergers and from exotic sources, such as cosmic string interactions 
(jDamour fc Vilenkin 20051 : IPshirkov fc Tuntsovll2010ft . Precision timing has also been proposed as a means 
for dete cting low-mass primo rdial black holes (< 10 -3 Earth masses) that make close flybys of the solar 
system (|Seto fc Cooravl 12007ft . Individual GW sources may be detectable if they are sufficiently massive 
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and nearby ()Finn fc Lommer] 12010). No matter the source, the targeted signals are very weak, requiring 
timing precisions along with stability of the pulsar clocks < 100 ns for a large number of MSPs over many 
years. 

We present a measurement model that allows an end-to-end assessment of timing precision. We distinguish 
between phenomena that affect the time-tagging of pulses — TOA estimation — and the astrophysics of 
the pulsar clock. An illustration of this distinction consists of pulses that can be time-tagged to arbitrary 
precision, such as delta-function pulses with infinite signal-to-noise ratio, but where stochasticity in the 
pulsar spin rate and in interstellar delays introduces a random component into the arrival times. We 
thus separate terms in the timing equation that result from the quality of the astrophysical clock and 
over which we have no control (except by choice of pulsars) from those that can be mitigated through 
choices of observational parameters and by appropriate post processing. These include phase jitter intrinsic 
to the pulsar, instrumental effects, and all the strongly chromatic interstellar propagation effects. We 
acknowledge but do not analyze timing errors from instrumental polarization, time transfer, and solar 
system ephemerides, which are discussed in the literature and likely will be improved to the point of 
becoming secondary factors in timing precision. 

Timin g perturbations that are interstellar in origin are one of the main topics of this paper. [Armstrong] 
(|1984l ) first pointed out that stoc hastic dispersion-me asure variations could limit timing precision. A timing 
and interstellar scattering study ([Cordes et al.lll990[ ) of the MSP B1937+21 (J1939+2134), indicated that 
diffractive interstellar scintilltions (DISS) cause fast timing perturbations and that epoch dependent inter- 
st ellar variations also occur. Diffractive TOA variations were also identified in observations and simulations 
bv lColes et ail ([2010D . iFoster fc Corded (|l990h assessed the roles of dispersion measure and refraction-angle 
variations and concluded that non-dispersive interstel lar effects need to be mitigated i n precision timing 
programs. This was followed by a similar study by iHu et al.l (jT991i) . More recently, lYou et ail (|2007|) 
demonstrated the importance of correcting for dispersion measure variations in precision pulsar timing. 

Our analysis leads to discussion of possible methods for removing most of the interstellar effects along with 
some of those intrinsic to the pulsar that affect TOA estimation. If measurements were made with the same 
telescope and the same intrumentation in a long timing program and that, hypothetically, the ISM caused 
chromatic perturbations that were epoch independent, much of the mitigation effort would be unnecessary. 
However, ISM effects are inherently stochastic and PTAs are expected to make use of multiple telescopes 
with instrumentation operating in different frequency bands that likely will change with time. Consequently 
mitigation of ISM effects must be done carefully. 

Our studies also apply to cases where interstellar scattering limits the timing precision to milliseconds or 
worse, which is relevant to any pulsars in the Galactic center even when observed at frequencies above 
10 GHz. 

In a separate paper ([Shannon fc Cordes Il2010l hereafter Paper I) we assess the role of spin noise in preci- 
sion pulsar timing, which we mention here only briefly. A companion paper (Shannon & Cordes 2010, in 
preparation; hereafter Paper III) presents simulations that corroborate and probe more deeply the inter- 
stellar plasma effects we discuss here analytically. We defer to another paper a consolidated analysis of the 
requirements of a PTA for making detections of GWs. 

In § 2 we summarize prominent physical contributions to arrival times. In § 3 we focus on achromatic and 
weakly chromatic timing perturbations. § 4 concentrates on strongly chromatic effects. In § 5 we consolidate 
various physical terms into a signal model for use in fitting arrival times. In § 6 we discuss the efficacy of 
different schemes for correcting arrival times for chromatic interstellar contributions. We summarize our 
conclusions in § 7. Appendices [A] and [B] present details of several contributions to arrival-time estimation 
errors. Appendix [Cl considers the apparent frequency dependence of dispersion measures that results from 
scattering. Appendix IDl gives the details of a least-squares fitting formalism. 

2. TERMS IN A PHYSICAL MODEL FOR PULSE PHASE 

Timing equations have been presented (e.g. iBacker fc Hellingslll986t lEdwards et al.|[2006l ) that give a com- 
prehensive inventory of deterministic contributions to the overall arrival time and, indeed, motivate timing 
studies in the first place. However, they have only incompletely addressed the timing error budget. Im- 
provements in timing precision are always being strived for that allow new effects to be measured. We 
describe salient terms in the timing budget that must be considered in any precision timing analysis. We 
first categorize them in terms of their physical origin. Then we consolidate terms according to their statisti- 
cal properties and wavelength dependence. Table [T] lists a comprehensive (but by no means exhaustive) set 
of terms and defines notation for deterministic and stochastic components. We designate the chromaticity 
of each effect, the spectral signature, and whether there is a correlation between different pulsar lines of 
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sight. Entries in Table [T] are labelled as astrophysical ('A' in column 2) if they determine the arrival time. 
Others are designated with 'T' if they contribute to TOA estimation errors. Some effects of course fall in 
both categories. Here we describe briefly some of the systematic astrophysical factors that determine the 
mean arrival time. A detailed discussion of interstellar perturbations is given in later sections. 

Pulse phase is dominated by the spin rate of the NS and its deceleration by magnetospheric torques. The 
stochastic part includes glitches and timing noise, the latter ty pically showing low - freqen cy components with 
a red power spectrum and more rarely a bandpass spectrum (|Hobbs et a,l.ll2009l I2010T) . Timing noise may 
comprise torque pulses or step functions that arise f rom effects internal to the neutron star (superfiuid-crust 
interaction) or within the magnetosphere. Recently iLvne et al. I (|2010t ) have shown that torque fluctuations 
in the form of jumps between discrete states at quasi-periodic intervals comprise much (but not all) of 
the timing noise in long-period pulsars. Spin fluctuations clearly are a nuisance for the use of pulsars as 
precision clocks, as is well known. The power spectrum for spin variations in some pulsars is degenerate 
or ne arly so with that expected from some sources of cosmological GW backgrounds (|Shannon fc Cordes I 
2010 ), so un l ess tim ing noise can be removed, it will play a significant role in the choice of pulsars for PTAs. 
Lvne et al. I (|2010f) show that distinct pulse shapes are correlated with the discrete torque states in some 
objects. The correlation must be close to unity to remove most of the timing noise and so far it is not clear 
that the correlation is large enough for timing noise to be ignorable in forecasts of timing precision. 

Pulses emitted from the rotating magnetosphere are proxies for spin phase of the NS (with due allowance 
for any doppler shift from NS motion), but emission regions in the magnetosphere are not synchronized 
perfectly with the NS. Arrival times are (usually) calculated by fitting a template to pulse profiles obtained 
by averaging N pulses. Individual pulses are well known to vary in amplitude and pulse phase by large 
amounts, but by mechanisms that are not well understood. However they likely i nvolve incomplete filling of 
the o pen-field line region in the magnetosphere by coherently radiating plasma (|Ruderman fc Sutherland! 
If 9751 ). Radio emission may occur at different altitudes in response to time variable pair plasma densities, 
which will change arrival times in accord with light-tr avel times, angular aber ration, gravitational bending 
of ray paths, and any refraction in the magnetosphere ()Barnard fc Aronslll986|) . Arrival time perturbations 
from the finite number of pulses combined with intrinsic phase jitter imply that there is a correlation of 
TOA error with pu lse shape. These changes in shape are much different in origin from those identified by 
ILvne et al. I (|2010f) and signify changes in pulse phase rather than torque. Nonetheless, a sufficiently high 
correlation again implies another opportunity for correcting TOAs, which we explore elsewhere. 

The evolution of pulse shapes with frequency includes changes in widths, separations and amplitudes of 
individual pulse compo nents that comprise av erage profiles. Pulse widths typically scale with frequency as 
W oc V ' x with x < 0.3 (jMitra fc Rankinll2002l) . so observations that strive for TOA precisions of 1% of the 
pulse width or better require mitigation of pulse shape changes. A key issue is the alignment of frequency- 
dependent profiles so that TOAs are consistent between frequencies. One prescription is to identify fiducial 
pulse phases at differen t frequencie s so that the resulting TOAs match those expected from the cold-plasma 
dispersion law, cx v~ 2 (|Craftl[f970t ). This procedure will fail at the highest timing precisions because there 
are other frequency-dependent, interstellar time delays that scale differently with frequency (|Romani et al.l 
I986; ICordes et al1ll986tlFoster fc Cordeslfl990[) . A recent discussion of the inter action betwe e n pul se shape 
evolution with frequency and estimates of dispersion measures may be found in lAhuia et al.l (|2007l) . While 
the focus of this paper is on frequency-dependent plasma effects, a complete analysis must include all pulse 
shape changes. 

Other entries in Table [T] correspond to deterministic contributions for orbital motion, gravitational lensing, 
and astrometric errors associated with referencing TOAs to the solar system barycenter (SSBC). Any error 
in th e location of the SSBC from uncertainties in planetary perturbations or from passing sub-stellar objects 
(e.g. ISeto fc Cooravl [20071 ) will produce a perturbation having a dipolar signature, which we designate as 
A^Nowt,ssBC- Stochastic perturbations from any orbital debris (R. Shannon e t al., in preparation) and 
instru mental effects from time transfer and polarization calibration errors (e.g. lBrittonll2000l ; Ivan Stratenl 
2009) are included in the table. 

While all listed phenomena must be confronted in precision timing programs, this paper focuses on those 
that affect estimation of individual TOAs, including pulse shape stability, additive noise, and especially 
delays that arise from dispersion, scattering and refraction in the interstellar plasma. 
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a A = astrophysical, T= timing estimation error 
b a = achromatic, C = strongly chromatic, c = weakly chromatic 
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d U = uncorrelated between different pulsar lines of sight, C = correlated 
c Includes clock errors and Earth spin variations 
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3. ACHROMATIC AND WEAKLY CHROMATIC TOA PERTURBATIONS 

The consolidated stochastic and achromatic term (including weakly chromatic contributions listed in Ta- 
blet!]) is 

At s , A = Aij + AUst + Ai s/N + Atpoi. (1) 

In the following we are most interested in terms that appear even when astrometric contributions and 
instrumental polarization are dealt with perfectly. 

Assuming that astrometric and polarimetric errors are negligible after appropriate processing, a simplified 
stochastic, achromatic term that includes only pulse phase jitter and radiometer noise, is 

At SA = Atj + Ai s/N . (2) 

Both terms affect the shape of the pulse profile formed from a finite number of pulses and thus contribute 
to the random error in the arrival time. 

3.1. Radiometer Noise 

Radiometer noise adds to the p ulse shape measured in a ny of the Stokes parameters. The resulting TOA 
error standard template fitting ()Downs fc Reichlevlll983f) is approximately VF/SNR, where W is the pulse 
width and SNR is the signal to noise ratio of the pulse. The frequency dependence of the sky background is 
strong (T s k y oc v~ 2 ' 7 ) at frequencies below about 0.5 GHz. However, for higher frequencies usually used for 
precision timing, the frequency dependence is weaker. The signal-to-noise ratio involves the pulsar's radio 
spectrum that typically also declines with frequency, canceling much of the variation of the sky brightness 
temperature with frequency. We designate radiometer noise as weakly chromatic in Table [1] an d bundle 
it wi th stochastic achromatic terms. In some cases, however, the pulsar's spectrum is flat (Lori mer et al.l 
1995), so the signal-to-noise ratio can be strongly chromatic. 

In Appendix [X] we derive the minimum TOA error that applies when the conditions for matched filtering 
are satisfied, namely that the measured pulse shape comprises an invariant templae shape added to noise; 
our analysis assumes that the noise is white. For a sum of N = 10 6 A6 pulses, is 

At s/N = 1 ^sW^iV-^SNR^A/W) 1 / 2 = 0.71 fis ( V /1A GHz)- Qp P^W^B N^ 1 / 2 (S sys /S 140 o) -(3) 

In the first form, SNRi is the signal-to-noise ratio of a single pulse (defined as the ratio of pulse peak 
amplitude to rms noise) and W is the effective pulse width defined in Appendix [XI expressed here with 
a fiducial value of 1 ms. For a Gaussian pulse shape the effective width is W = Wfwhm / V27r In 2 « 
0.69Wfwhm- The quantity A is approximately equal to the time resolution. In the second form, which 
shows explicitly that Atg/N is independent of A, S sys is the system temperature expressed in Jy; S1400 and 
a p are the period-averaged flux density in mJy and the spectral index, respectively; the radio frequency v 
and bandwidth B are in GHz; and the spin period P is in ms. 

The appropriate values for pulse width W and SNRi include any distortions that broaden the pulse and 
reduce its amplitude. Residual dispersion smearing is one such effect, although it is negligible if pulses are 
coherently dedispersed. Broadening from interstellar scattering of distant pulsars observed at low frequencies 
makes the pulse asymmetric and thus non-Gaussian, so the expression is only approximately correct. 

Another way that the ISM affects the timing error is through DISS of low-DM pulsars observed at high 
frequencies, which can cause the apparent flux density to be much smaller or much l arger than average. 
TOAs are sometimes obtained only during scintillation maxima (e.g. iHotan et al.ll2006l) which can enhance 
the SNR by two orders of magnitude, reducing the noise-induced error below other contributions. 

3.2. Pulse Jitter and Amplitude Modulations in Known Pulsars 

Al l well-studied canonical p ulsars (CPs) show phase jitter and amplitude modulations like that in Figure 4 
of [McLaughlin et al. (2002). Individual pulse amplitudes vary by 100% (or more) and pulse phases jump 
by of order the single-pulse width. Empirically, the average pulse shape of a pulsar is as much deter mined 
by the PDF of the phase jitter as it is by the shapes of individual pulses, a description first made by ICraftl 
(|1970[ ). The commonly- used intensity modulation index to/ (the rms intensity divided by the mean as 
a function of pulse phase), is typically of order unity and quantifies both the amplitude modulation and 
the phase jitter. Phase jitter appears to be statistically independent between pulses in some objects, while 
others show systematic drifts over several to many pulse periods, yielding an effective number of independent 
pulses Ni < N in an iV-pulse average. We distinguish jitter from sustained "mode" changes in pulse shape 
displayed by some pulsars, where pulse shapes switch between two or three preferred shapes. We quantify 
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phase jitter alone; if mode changes are unaccounted for in TOA estimates, they will only exacerbate the 
timing estimation error. We use the dimensionless quantity fj to quantify the rms phase jitter in terms of 
the intrinsic pulse width, Wi. 

Objects other than CPs also show phase jitter. The Crab pulsar shows giant pulses with amplitudes several 
hundred times larger than the m ean pulse amplitude; giant pulses are int r inically narrow (1 — 10 fis) 
but show TOA jitter ~ 100 us (|Lundgren et all 119951; ISallmen et al lll999t ICordes et all 12004) . Giant 
pulses from the MSP B1937+21 show fj « 0.2 - 0.5 (jKinkhabwala fe Thorsettl 120001) . Howe ver, in stark 
contrast, the ordinary (non-giant) pulses from this pulsar apparently show negligible jitter (jJenet et all 
120011; I Jenet fc Gilll2004). Other MSPs show jitter similar to t hat of CPs, including J0437-4715, the brightest 



known MSP (j Jenet et aTlll998l) . lEdwards fc Stappersl (|2003[ ) show that the autocorrelation function (ACF) 
of the average profile for the recycled 41-ms pulsar J1518+4904 is about twice the width of the ACF of 
single pulses (their figure 9), an effect that can be accounted for only with phase jitter. 

Jitter is related to pulse-shape stability that has been qua ntified using the cross correlation coefficient p(N ) 
between a template profile and sub-averages of N pulses (jHelfand et allll975t iRathnasree fc Rankinlll995l) . 
These analyses typically show 1 — p declining as TV -1 for some ranges of N (as expected for statistically 
independent phase jitter between pulses) while in some objects the trend is flatter for small N, signifying 
correlations between pulses. It can be shown that the correlation coefficient is 

P (N = !)=[(! + mj)(Wu/W a )] 1/2 = (l - fj) 1/4 / (l + mf ) 1/2 , (4) 
where Wjj is the characteristic width of the template profile, W a is the characteristic width of an individual 
pulse and we have used the relation, fj = l — W^/W^j. Typically with m a = 1 and Wjj — 2W a , we expect 
p(N = 1) = 1/2. For large N we have p(N) « 1 - (l/2iV)(l + ™>l)(Wu /Wa)- Correlation values reported 
in previous stability analyses are consistent with these estimates, indicating that fj ~ 1/3 to 1/2 for most 
pulsars. 

3.3. Timing Error from Pulse-phase Jitter 

When TOAs are calculated from averages of N pulses, any given profile will depart stochastically from a 
long-term average, even when the SNR is large. Template fitting therefore yields a jitter-induced error 
that adds quadratically to the minimum expected from matched filtering . The role of jitter in pulse-shape 
stabi lity and timing has been discussed in the past (|Cordes fc Downslll985tlD'Alessandro et al.lll993HCordesl 
119931 ). but little overall attention has been given to its role in precision timing. 

The TOA error cx Wi/y/W is pulsar dependent. (Appendix |A")) . Because single pulses from many pulsars 
appear to be broadband with modest frequency dependence, the jitter timing perturbation in a typical data 
set is only slightly dependent on frequency and can be viewed as a stochastic element of the systematic 
profile evolution in frequency described earlier. For simplicity, in some of our analysis we will consider the 
jitter effect to be stochastic in time but constant in frequency. For Gaussian-shaped pulses and a Gaussian 
probability density function (PDF) for pulse phase jitter, 

At, = 0.28^ W itms N^ 2 ( 4) ( ■ (5) 



,1/3, 

Figure [1] shows the rms TOA error caused by different levels of phase jitter from simulations. The left-hand 
panel shows a sequence of pulses that displays pulse-phase jitter with fj = 1/3. The right-hand panel 
shows the rms arrival-time error for pulse profiles formed from a fixed number of single pulses (10 3 ) as a 
function of signal to noise ratio and for different values of fj. For fj = 1/3, the rms error is double that 
from additive noise once the single-pulse SNR exceeds a few tenths, which is consistent with a comparison 
of Eq. O [5] 

In Figure[2]we show a jitter analysis for the MSP J1713+0747. Data were obtained at the Arecibo Observa- 
tory at 1.38 GHz using a baseband recording system, coherently dedispersed to yi eld high-time-res olution, 
and averaged synchronously using a timing model from P. Demorest and D. Nice (|Demo rcst 20071, private 
communication). Details of the data acquisition and other analyses will be described elsewhere (Shannon 
et al., in preparation). We compare TOA errors expected from radiometer noise alone to those actually 
measured. First we take a high S/N template profile, add white noise, and then calculate the arrival time. 
Repeating over a large number of realizations and for different amounts of noise yields the open circles 
shown in Figure [21 We use Eq. [3] and calculate the effective width using the template profile and Eq. IA2I 
The numerical rms TOA error agrees very well with the predicted error. Actual data show excess TOA 
variations. We calculated TOAs for an initial set of 10-s averages (2188 pulse periods) and then for 20-s, 
40-s averages, etc. A jitter parameter, estimated assuming to/ = 1, 

fj = [N (4oa(N) - a TF (N) 2 ) I (1 + mf)] 1/2 « 0.4 (6) 
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fj = 0.333 f G = 0.000 S/N = 100.000 FWHM = 0.050 
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Single Pulse S/N 



Fig. 1. — (Left) Simulated pulse sequence of showing pulse phase jitter with fj = 1/3 and a unity intensity modulation 
index. The net duty cycle is 0.05. (Right) Arrival time error vs signal to noise ratio (SNR) for simulated data with different 
amounts of pulse phase jitter, as labelled. TOAs were calculated from averages of 10 3 pulses to yield a SNR for each profile 
that is 10 3 / 2 larger than the single-pulse SNR. Each plotted point is based on 500 realizations. 



accounts for the excess timing error over that expected from template fitting (TF) to a profile consisting of 
the template combined with additive noise to yield the same SNR as the real data. 




0.01 \= I I I L_d 

10 10 2 

S/N 

Fig. 2. — Arrival time analysis for the MSP J1713+0747 at 1.4 GHz. The points are calculated for pulse profiles calculated 
from different numbers of pulses (N, top axis) yielding increasing signal-to- noise ratios (bottom axis). The lower set of points 
results from calculating TOAs from simulated data formed by adding a high SNR template to gaussian noise. The upper 
set of points is from fits to actual data obtained with different integration times. The short-dashed curve is the fit to the 
points. The long-dashed curve is the prediction using Eq. IAUA2I which assumes that data comprise identical pulse shapes 
with additive noise. The two curves are actually coincident but have been displaced by ±0.01 in the log for clarity. The 
actual data have higher TOA errors than simulated values because pulse-phase jitter causes pulse shapes to differ from the 
template shape. A jitter parameter fj fa 0.4 accounts for the difference. 
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3.4. Jitter vs Noise Dominated Timing Errors 

We compare jitter an d noise contributions to the TO A error for the known sample of pulsars by using 
cataloged parameters (M anchester et al~ll2005f) and assuming reasonable values for amplitude modulation 
(mi = 1) and phase jitter (fj = 1/3). Figure |3] gives the fraction of MSPs for which phase jitter dominates 
the radiometer noise error. We assume that jitter is statistically independent between pulses so that Ni = N, 
which may underestimate the TOA error. Fractions are calculated by counting objects for which Atj > 
A^s/n for different system temperatures (in flux-density units, Janskys) that correspond to 100-m class 
telescopes (e.g. the Effelsberg and Green Bank Telescopes, the Expanded VLA, and the Parkes and Sardinia 
telescopes) , Arecibo-class telescopes ( Arecibo itself, the Chinese FAST telescope, and Phase 1 of the Square 
Kilometer Array, SKA), and Phase 2 of the SKA, which is projected to have S sys of order ten times smaller 
than for Phase 1. Results are shown for frequencies from 0.3 to 5 GHz using either the cataloged spectral 
index of the pulsar's flux density or using a power law oc is~ 2 . 

The results indicate that usage of more sensitive telescopes (whether by larger collecting area, lower system 
noise, or greater bandwidths) will lead to timing measurements that are limited by phase jitter. To increase 
the timing precision, the only recourse is to make use of longer integration times per TOA. It may also be 
possible to reduce the timing error from phase jitter if the resulting pulse shape variations are correlated 
with the timing error. A method based on principal component analysis will be reported separately (Cordes 
et al., in preparation). 
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Fig. 3. — Fraction of known pulsars with jitter timing error larger than that from radiometer noise, i.e. Afj > Afg/ N . 
The curves are labeled by the system equivalent flux density, S sys , at 1.4 GHz. S ByB varies with frequency by adding the 
Galactic synchrotron contribution as 20 K (u/0A08 GHz) -2,75 as quantified by Cortes (2007). Calculations were done with 
10% bandwidths below 1 GHz and bandwidths of 0.4 and 0.8 GHz at 1.4 and 2 GHz, and 1 G Hz bandwidths at 3 GHz a nd 
higher. Each MSP with P < 20 ms that had tabulated flux densities in the pulsar catalog of [Manchester ct al. (2003) is 
included. Given that only 82 MSPs in the current catalog satisfy the criteria, a cubic fit was done to the computed fraction 
to plot each curve. 




4. STRONGLY CHROMATIC TOA PERTURBATIONS 

Stochastic and strongly chromatic timing perturbations include contributions from plasma dispersion, Fara- 
day rotation and interstellar scattering (ISS), 

At s , c = At DM + At RM + Atiss- (7) 

The dispersion term receives contributions from the ionosphere, interplanetary medium (IPM), and the 
ISM. Any extragalactic pulsars will also involve the intergalactic medium and the ISM of their host galaxy. 
The ISS term comprises six separate terms 

Atiss = AtpBF + A<AOA + At A OA,SSBC + A^pbf + Ai 5 pBF,RISS + Ai DMil/ . (8) 
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These are respectively the delay associated with diffractive pulse broadening, the delay from angle-of- 
arrival (AOA) variations, the error in referencing TOAs to the SSBC due to AOA variations, the error from 
fast variations in the pulse-broadening function (PBF) that describes diffractive pulse broadening, slow 
variations in the PBF induced by refraction, and the perturbation of DM that results from averaging over 
the scattering cone. We discuss these terms in detail in § 14.31 

The relative amplitudes of ISS effects depend on the spatial properties of the ISM. The standard approach 
is to characterize small-scale variations in electron density (< 1000 AU) — those that diffract and refract 
radiation from point sources — with a wavenumber spectrum and to consider large-scale variations (e.g. 
cloud-like, sheet-like and filamentary structures) as modulations of the spectrum. Small scale variations 
largely appear to be consistent with a Kolmo gorov spectrum and in many cases the variations involve 
anisotropic irregularities with large axial ratios ( Cordes et al J 120061 ; iBrisken et al.ll2010|) . Regardless of the 
detailed shape of the wavenumber spectrum , so long as it is broad, three regimes of scattering can be 
defined (Rickctt 1990; Cor des fc Lazio| [l991). When the phase perturbation by the scattering medium is 
much smaller than a radian across one Fresnel scale ~ \/\D c g , where A is the wavelength and D c g is a 
characteristic distance, a point-source will appear largely as a point source with a small halo around it. 
When the phase perturbation is large, the scattering is strong. All of the source flux will be scattered and 
the scintillation bandwidth Ai>d <C v. So long as the scattering length l& = X/9d is larger than any minimum 
(inner) scale li in the wavenumber spectrum, the scaling of Avd with v is a function of the spectrum's slope 
(if shallow enough). For super-strong scattering, Id < t\ and the scaling laws generally will be different 
than in the strong-scattering regime. In particular, they will conform to a medium with a square-law phase 
structure function. Electron density fluctuations appear to exist on a wide ran ge of scales from ~ 10 3 km 
to Galactic scales so, accordingly, there is a large range of expected time scales ([Armstrong et al.lll995l) . 

When scattering is weak, only a fraction of the pulsar's flux is scattered, the remainder appearing as a point 
source. In this case, the delay of the pulse is negligible because even the scattered portion is delayed by only 
(2-7ri/) -1 R3 0.1 ns. The transition from weak to strong scattering occurs approximately when Avd/v = 1 or 
^trans = C\/ (2n~fd)- Using Td = TdQV~~ AA , the broadening time at 1 GHz and with v in GHz, (the Kolmogorov 
scaling when there is neglibile inner scale), we obtain ^t ians = (2nTdo/Ci) ' = 12.1 GHz T d0 (^s) a29 . From 
Figure 21 which shows values of Td at 1 GHz plotted against DM, the transition frequency is about 2.3, 12 
and 24 GHz for pulsars with DMs of 10, 50 and 100 pc cm" 3 . 



4.1. Dispersion Measure Variations 
For a total DM from all contributions along the line of sight (LOS), the group delay is 

*dm - odm DM v- 2 (9) 

where odm = 4.15 ms for DM in standard units of pc cm -3 . Fluctuations JDM(t) affect the arrival time itself 
and, if large enough and unaccounted for, can distort the pulse shape and affect the TOA. In the following we 
ignore such distortion because it is addressible in the data processing. DM variations produced by structure 
in the ISM are stochastic in time but deterministic in frequency (however, see § I4.4[) . Measured variations 
are seen on a variety of time scales that are consistent with a broad spectrum of electron-den s ity variations 
in th e ISM (1 Armstrong 111984 ICordes et al.Ml990t iRickettl [l990t iPhillips fc Wolszczad Il99ll: iBacker et all 
1993: lKaspi et al.lll994t ICognard fc Lestradelll997|) . Amplitudes depend on the distance and velocity of the 



pulsar but maximum rates of change ~ 10 -3 pc cm" 3 per year. The te mporal fluctuation spectrum for DM 
has the same form as the electron-density spectrum ([Armstrong 1119841) and has large fluctuations on longer 
time scales. During a typical observing session at a given epoch (e.g. < 1 hr), DM variations are too small 
to be detected with known pul sars. As is well known, however, precision timing requires removal of DM 
variations between epochs (e.g. lYou et al. I [20071) . 

The precision needed for DM is stringent. The simplest case is when the only other TOA perturbation 
is from white-noise measurement error. We write the TOA as the sum of achromatic, dispersive, and 
white-noise terms: 

t v = too +*dm(i>) +*w(^), (10) 

where is the achromatic arrival time at infinite frequency (which would include spin noise and gravi- 
tational wave contributions). Here we ignore any frequency dependent profile "evolution" intrinsic to the 
pulsar. For a two-parameter (too and DM) least-squares fit across a frequency range [v2, vi] sampled in 
n v channels, the error on the achromatic TOA is — S4/(sqS4 — s§), where s p is a weighted sum over 
frequencies defined in Eq. ID1I For equal errors in all channels er too depends only on the frequency ratio 
r = v\jvi > 1, as expected. An alternative method would use a value for DM obtained from ancillary, 
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non-contemporaneous measurements to correct TOAs rather than fitting for DM from the TOAs. For this 
procedure to be at least as good as fitting for DM, the error in <5DM must satisfy 

v z 

<5DM < l(T 3 - 38 pc cm" 3 -i(l - B/^i) 3/2 ct b (» (11) 

B 

for v\ and B = V\ — v-i in GHz and a fiducial radiometer noise TOA error ub = 1 fxs (c.f. Eq. [3]). In 
this expression, we assume that errors are independent of frequency and emphasize that this is just the 
radiometer noise contribution to the error. 

Measured DM changes are at levels expected for interstellar electron-density variations that have a power- 
law wavenumber spectrum. For a power-law spectrum C 2 X (wavenumber) - ^, where C 2 is simply a spectral 
coefficient in the wavenu mber spectrum, the structure function f or DM has an ensemble average, denoted 
by angular brackets (e.g. lCordes fc Ricket tl 1 1 9 9 8t lYou et al.ll2007l ). 

Ddm(t) = {[DM(t) - DM(t + r)] 2 ) = f a J ds |V effjL (s)rC 2 (s) « ^fjSM (rV p ) Q , (12) 

where a — j3— 2, f a is a constant, V c ff^ is an effective velocity, and SM, the scattering measure, is the LOS 
integral of C 2 . The approximate equality holds for constant C 2 and when the pulsar transverse velocity 
dominates the effective velocity V e g ± . For a Kolmogorov medium, / s / 3 = 88.3, and using standard units for 
scattering measure (kpc m~ 20 / 3 ), the rms difference in DM between two epochs separated by r (in years) 
is 

<tdm(t) = [(l/2) J D DM (r)] 1/2 « ltr 3 - 47 pc cm" 3 (sM/l(T 3 - 5 kpc m" 20 / 3 ) ^ (r yr V P100 ) 5/6 (13) 

for pulsar velocity in units of 100 km s _1 . This rms is comparable to measured changes in DM for nearby 
pulsars that have SM sa 10 -3 5 kpc m -20 / 3 . From Eg. [TBI and [TT1 if we require a timing error less than 0.1 fis 
for 3ctdm(t)) then the DM and TOA measurements need to be made within r < 8 days. This result is in 
accord with those determined numerically in Paper III. 

4.2. Profile Splitting from Faraday Rotation 

The TOA perturbation from a magnetized plasma includes a birefringent term that is equal and opposite 
for the two hands of circular polarization, £rm = ±0.18 ns i/~ 3 RM for RM in standard units of rad m~ 2 . 
The empirical manifestation is profile splitting of the two polarizations. For the extreme case of high-RM 
pulsars observed at low frequencies, the profile splitting produces a pseudo-circular polarization that would 
cause a sign-changing Stokes parameter V = II — Ir even when the circular polarization intrinsic to the 
pulsar is zero. If uncorrected, the total intensity profile I = Il+Ir is broadened slightly but symmetrically 
relative to the true TOA, so the TOA error is not biased one way or another. Correction of this effect, if ever 
necessary, needs to be done prior to summing the two hands of polarization and can be incorporated into 
the coherent dedispcrsion process. Values of RM for Galactic lines of sight can range up to 10 3 rad m -2 but 
for objects within 1 kpc, typical RMs are a few tens of rad m -2 or less. The RM perturbation is therefore 
negligible compared to a 100-ns timing error unless distant pulsars with large RMs are observed at low 
frequencies. 

4.3. Interstellar Scattering and Refraction 

The interstellar timing perturbations in Eq. [5] result from refraction and multipath propagation caused by 
scattering from small-scale electron-density variations. The resulting timing perturbations have distinct 
time-frequency characteristics: 

1. Variations in dispersion measure from frequency-dependent spatial averaging of the phase 
perturbation (AtDM,i>)i 

2. The mean time delay associated with multipath broadening of pulses from diffraction (A£pbf); 

3. Fast, stochastic changes in the diffractivc delay associated with the finite number of "scintles" in 
the sampled time-frequency plane (A^pbf); 

4. A time delay associated with the stochastic refraction angle of arrival (AOA, A^aoa); 

5. An AOA induced error associated with transforming measured TOAs to the SSBC (A£aoa,ssbc), 
which arises from the difference between the assumed and actual direction to the pulsar; and 

6. Modulations of the mean time delay from wavefront curvature induced by large-scale variations, 
which also causes refractive interstellar scintillation (RISS; A£,5pbf,riss)- 
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The first effect is the modification of DM by the frequency-dependent spatial averaging of the electron- 
density along different propagation paths. The second (typically dominant) effect is the delay introduced 
by the pulse broadening function (PBF) that is convolved with the true pulse shape. The PBF is causal, so 
the delay is always positive. It is also strongly frequency dependent. The PBF is stochastic from refractive 
variations and from the scintillation fluctuations contained in the finite bandwidth and integration time 
used to form a pulse profile. The angle of arrival (AOA) of the centroid of the bundle of multipath 
"rays" varies from refraction in the ISM and adds another frequency-dependent delay with a different 
freque ncy dependence. Final ly, the AOA induces an error in the correction of topocentric TOAs to the 
SSBC (|Foster fc CordesllT990l) . 

We characterize the second and last three effects as a mean delay that changes slowly with time (hours to 
years). The third effect differs from the others because it involves interstellar diffraction features that are 
statistically independent between data sets used to calculate individual TOAs. We give a detailed analysis 
in the next section. 

For a Kolmogorov medium in the strong-scattering regime, refraction angles are much smaller than diffrac- 
tion angles, so the TOA perturbations included in Eq. [8] may be ordered as 

AtDM.i/ ^ AiAOA < AtAOA,SSBC ^ At,5PBF,RISS ^ Ai^pBF < A£ PB f (14) 

because the two smallest terms depend quadratically and linearly on the refraction angle, 8 r , while the other 
terms depend on powers of the diffraction angle, 9d- Scaling laws and estimates are given in § 14.31 Alternative 
orderings of terms may apply for media with enhanced refraction from non-Kolmogorov structures. 

The scaling of the pulse broadening time with frequency indeed shows departures for some objects from the 
scaling expected for a Kolmogorov medium (Lohmer et al. 2003). These departures either indicate actual 
departures from a Kolmogorov wavenumber s pectrum or the breakd ow n of the assum p tion o f an unbounded 
screen, which can give anomalous scalings ()Cordes fc Laziol [2001[ ). lLohmer et al.l (|2004f ) conclude that 
anomalous scattering is associated with distant objects whereas nearby objects generally sho w Kolmogorov 
scalin gs (or nearly so) . One such scaling is for the scintillation ban dwidth A.Ud with frequency (jCordes et al.l 
1985). Another is the slope of the structure function for DM (e.g. lYou et al.ll2007f ). 

4.4. Frequency- dependent DM Variations From Multipath Screen Averaging 

Consider a simple screen model for the diffraction and refraction. The diameter of the scattering cone at a 
phase screen located a distance D' = D — D s from Earth for a pulsar distance D is t s « ID'Od- Measurement 
of DM toward a source therefore involves spatial averaging over an area ~ t 2 s . This causes the measured 
DM to differ from that for the direct line of sight through the screen. Since 9d oc v~ 2 2 or so, this implies 
that the difference in DM is strongly frequency dependent. By defining a weighting function W(x, v) that 
is defined by the scattered image expressed with a spatial argument (i.e. x = D'9) and thus has a width 
~ £ s , it can be shown that the rms difference in DM is a weighted integral over the phase structure function 
and 

SBM re 10~ 105 pc cuT 3 SMD' 5/ V 11/6 (15) 

for SM in standard units, D 1 in kpc, and v in GHz. The resulting TOA perturbation, which is stochastic 
with characteristic time scale £ S /V e s ± (the same time scale as for RISS), has rms 

A W « 0.12 ,s D*'%-^ (-^ m ) • (16) 

While small, this effect results from the dispersion measure effectively being a function of frequency. It 
is a steep, declining function of frequency because the area of the scattering cone decreases rapidly with 
frequency. The frequency dependence is nearly degenerate with some of the other timing errors. 

4.5. Pulse Broadening Function 

Pulse broadening from multipath scattering in the ISM is described by a pulse-broadening function (PBF) 
that is convolved with the emitted pulse shape to produce the measured pulse. 

The ensemble average pulse-broadening function (PBF) has a characteristic width ~ D8 2 ,/2c where dd 
is the diffraction angle and D is an effective distance. The timing perturbation is proportional to the 
width of the PBF, which is strongly wavelength dependent (oc A ). The pulse broadening times 
demon strate that they are a strong function of both DM and frequency. The empirical fit from 
(|2004f) is logT d = -6.46 + 0.1541ogDM + 1.07(logDM) 2 - (3.86 ± 0.16) log v. The scatter about this fit is 
±0.65 in logr^. 
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We write the PBF as Pd(t) = Pdif) + 5Pd{t)- where p d is the ensemble- average shap43 expected from a 
particular medium and consider it to be a slow function of epoch. By constrast, dpd encapsulates fast 
departures from the ensemble average. Defining pd(t) to have unit area, the characteristic broadening time 
is 

T d = J dttp d (t) = r d + Sr d , (17) 

and its ensemble average as r d . 




1 10 100 1000 0.1 1 10 



DM (pc cm- 3 ) v (GHz) 

Fig. 4. — (Left) Pulse broadening time plotted against DM for pulsar data used to construct the NE2001 model 
l lCordes fc L azio 2002). The solid curve shows the parabolic fit given in the text for logr^. and the dashed curves show 
il.Oo" deviations from the fit. Open circles are inferred values using scintillation bandwidth measurements; filled circles are 
direct measurements of r^. Downward arrows denote upper limits. All measurements have been scaled to 1 GHz assuming 
Td <x i/— 4,4 and the right-hand scale for 2 GHz assumes the same scaling. The strong scattering regime is assumed, which 
requires ^> Ci/2nu or t,j > 2 X 10~ 4 £ts at 1 GHz. At lower frequencies, the entire range of in the plot satisfies the 
requirement. Leftward arrows designate maximum values of DM for which < 1 ms at 0.4, 1 and 2 GHz using the nominal 
parabolic fit. (Right) Contours of constant pulse broadening time, r^, using the fit shown in Figure [4] including the ±1(7 
uncertainty in the fit, which results in the widths of the shaded regions. A t<j oc v~ 4,4 scaling with frequency is used. This 
scaling is steeper than appropriate for some high DM objects but is typical for low-DM pulsars. 



The scintillation bandwidth and the mean pulse broadening time Td are related bjJ3 

2nT d Av d = Cu (18) 

where C± depends on all the properties of the medium (|Cordes fe Ricket"tlll998l ; lLambert fe Rickett | [2000). 
such as the wavenumb er spectrum, its thickness and location along the (LOS) and on its transverse extent 
(|Cordes fe Laziol 12001). For simple media, values of C± can be calculated. For a thin screen unbounded 
transverse to the LOS and with a circular Gaussian angular scattering function, the PBF is a one-sided 
exponential function with l/e time scale Td- In this case the mean time delay is also T d so that C\ = 
1. For media that are plausibly relevant to pul sar scattering, C± can vary by a factor of nearly two 
(jCordes fc Rickettl 1 19981 lLambert fc Rickett 111999ft . 



4.6. The Mean-Shift Regime 

i 

This might also be called a "snapshot" average using the language of Narayan & Goodman (1989), which takes into account 

that a statistically precise measurement at a given epoch is still not a good ensemble average because there are long-term 

refraction effects. Our approach is to treat refraction effects as a separate modulation imposed on a statistically precise average 

at a given epoch. 
2 

Our definitions of Au d and T d follow [Cordes & Rickett (199 j), i.e. Au d is the HWH M of the intensity ACF vs. frequency 
lag and is the mean delay of the PBF. This differs from Lambert & Rickett (1999), who use alternative definitions of t<j 
in a similar expression and thus quote different values of C\ than we use here. 
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A case of interest for nearby pulsars or observations at high frequencies is where the pulse broadening time 
is small compared to the pulse width, Td <C W but where the strong scintillation regime applies. This 
requires that the scintillation bandwidth Ai>d <C v or, equivalently, r d v 3> 1. We specify the mean-shift 
regime as 

C 

- — < r d < W or 0.2 ns v^L, «T d «lms W ms . (19) 

Zttv ' ' 

Under these conditions the TOA perturbation is simply the pulse broadening time calculated from Eq. [T7] 
using the instantaneous PBF (not the ensemble-average PBF). The TOA perturbation also depends on the 
number of scintles contributing to an observation. In many cases it is large and the TOA perturbation 
converges to some mean value. This requires that Av d <C B, where B is the bandwidth, or that Aiiss <C 
T, where T is the total time of the observation. Figure |4] yields the range of DM for which the shift 
approximation applies: < 1 ms for DM < 100, 250 and 400 pc cm~ 3 at 0.4, 1 and 2 GHz, respectively. 

One approach to timing precision is to require that Td be smaller than some specified rms timing error, <r max , 
which might be 100-ns or less for pulsar timing array applications but could be much larger for other timing 
purposes (such as pulsars in the Galactic center). If T d < cr max , pulse broadening might be ignored entirely. 
Alternatively, it might be corrected to a fractional precision e T , yielding a residual error e T Td < cr max . Using 
the empirical fit for T d vs. DM shown in Figure 2] we can then solve for DM(f) for different <7 max . Results 
are shown in Figure [5] The smallest pulse broadening effects occur in the lower-right portion of the diagram 
whereas objects with large DMs are strongly affected by pulse broadening. Many pulsars have steep spectra, 
so the choice of frequencies for a given object requires a compromise, as discussed in § [6] 

When the shift approximation holds, the TOA can be corrected by subtracting Td from nominal values 
(|Hemberger fc Stinebring| [2008). However, because by definition Yd is small, we cannot measure it directly 
at the frequency of the TOA measurements. There are several ways to estimate it, however. 



4.6.1. Correlation Approach 

The timing delay r d can be estimated by using Eg. 1181 where the characteristic scintillation bandwidth Avd is 
calculated as the half- width at half maximum of the intensity correlation function Y &i(bv) that is calculated 
from the dynamic spectrum I(i>, t). This approach requires knowledge of the shape of the PBF so the type 
of medium must be known to determine the appropriate value for C\. Alternatively, a characteristic time 
Tg 2 can be calculated fr om the secondary spec trum, the squared-magnitude of the Fourier transform of the 
dynamic spectrum (e.g. IStinebring et al.ll2001ft . which is related to the autocorrelation function of the PBF. 
The relationship of rs 2 to Td is also medium dependent and requires knowledge of their ratio that is analogous 
to knowing C\ . The correction of TOAs for the mean PBF delay is therefore problematic unless ancillary 
constraints on the scattering physics are known, namely whether a thin screen or extended-medium applies, 
whether the scattering region is bounded transverse to th e LOS, what the wavenumber spectrum is, etc. 
The s econdary spectrum provides much of this information ()Stinebring et a.L~ll2001t iHemberger fc Stinebring] 
12008ft . 

We identify the errors associated with this correction approach by considering Eq. [10] with an additional 
chromatic term, tc, added. For the case where tc consists solely of the pulse broadening term and the 
mean-shift regime applies, the TOA perturbation equals the pulse broadening time Td- The broadening 
time can be estimated from the scintillation bandwidth obtained at a frequency u', 

f d {v') = J* (20) 

where quantities requiring an estimate or an assumed value are hatted. The broadening time can be scaled 
to other frequencies using 

, / v \ ~A"pbf 

TdW) = W) {-) ■ (21) 

The estimation error of the chromatic term and thus on expands into terms dependent on individual 
errors, SC'i, SAvd, and <SApbf: 



$Td = Td - Td ~ T d 



SCi SAu d v 

— oApbf In — 

Ci Avd v 



(22) 



Based on estimates of scintillation bandwidths in the literature, the range of possible values of C\ , and the 
error in the exponent, <5Apbf, the combined error in Srd is unlikely to be less than 10% with this approach. 
However, a 10% error in removing pulse broadening implies that the shaded bands in Figure 2] will shift to 
the left by about a factor of 1.8 in frequency (0.26 in logi'). 
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4.6.2. Phase Retrieval Methods 

A second, empirical approach requires much less a priori knowledge about the scattering medium. It calcu- 
lates the complete PBF using only a measurement of the ACF of the PBF or, equivalently, the magnitude 
of the PBF's Fourier trans f orm. The ACF of the PBF can be obtained from the secondary spectrum, 
for example. Walker et al.l ()2008[ ) apply two-dimensional phase retrieval to the secondary spectrum and 
demonstrate the general principle. For the one-dimensional ACF with n T lags, there are n T functions that 
can produce the measured ACF. Many of these are acausal or include negative values and can therefore 
be ruled out. Those that remain include the function with minimum delay and its time reverse, which has 
maximum delay. Methods exist to calculate the minimum-delay PBF from the ACF, which is the unique 
member in the large family of PBFs con sistent with t he ACF that has its amplitude most concentrated 
toward the origin in a mean-square sense ([Scarglc 1981). While physically a minimum-delay solution is not 
demanded by scattering and refracting geometries, it may be the most probable result given that scattering 
and refraction angles have distributions peaked at zero. Also it may provide a more accurate correction 
compared to use of an idealized form for the PBF used to provide a value for C\ in the correlation method. 
Application of inversion techniques will be explored in a separate paper (Cordes & Shannon, in prepara- 
tion). Another approach (P. Demorest & M. Walker, private communication) investigates the phase of the 
wavefield directly to determine the PBF. 

4.7. Real World PBFs 

Time delays from pulse broadening presented so far represent ensemble-average results. TOAs, however, are 
determined from pulse profiles obtained as integrations over the time-frequency plane corresponding to the 
integration time T and bandwidth B. The T — B plane contains a finite number of independent fluctuations 
("scintles") of the diffractive interstellar scintillations (DISS), causing the instantanous PBF to differ from 
the ensemble avera ge shape, producing a statistical error in the estimated TOA, as first demonstrated by 
()Cordes et al.l 119901 ) . In addition, because individual TOAs are typically estimated from data sets with no 
overlap in T — B, the TOA error is statistically independent and thus has a white-noise spectrum. 

The TOA error arises from the summation of independent scintles in the T — B plane, the number of which 
can be large or small when a pulse profile is obtained. In the following we will use A^pbf to mean both 
the specific value applicable to a particular TOA and also to signify its rms value, the context making it 
clear as which applies. The rms PBF error is approximately 

AtiPBF = NvUaTd, (23) 
where Adiss is the effective number of "scintles." Previous work estimates -Adiss as 

A D iss = AW = (1 + TjtT/Atjss) (1 + V»B/Av d ) , (24) 
where Aiiss is the characteristic DISS time scale and Ai/^ is the characteristic bandwidth, as before. The 
filling factors f]ti^v are less than unity and are in the range of 0.1 to 0.3, depending on the definitions of 
the characteristic time scale and bandwidth. In the following we use approximate values for r\t and r\ v while 
keeping in mind that the factorization of Adiss is only an approximation. In another paper (Cordes & 
Shannon, in preparation) , we derive exact results for several cases that validate the form of Eq. 1231 

The scaling of the timing variation, A^pbf, with frequency is different for the four cases where T and/or B 
are larger or smaller than the characteristic scintle size, i.e. whether N t and/or N v are unity or much larger 
than unity and whether the scattering is strong or super-strong as defined in £14.31 In strong scattering 
T d oc z/ _i » where ( = 2/3/(/3 — 2) for /3 < 4. From Eq. [TH1 the DISS bandwidth therefore scales as ifi. 
For a Kolmogorov medium ((3 = 11/3), ( — 22/5. In super-strong scattering, the coherence length of the 
diffraction pattern is smaller than the inner scale and ( — > 4. Similarly, the characteristic DISS time scale 
scales as v K with k — /3/(/3 — 2) in strong scattering and k = 2 in super-strong scattering. The scaling 
depends on the type of medium (/?), whether the scattering is strong or super-strong, and on the bandwidth 
and integration time relative to the DISS bandwidth and time scale. For a Kolmogorov medium the index x 
in A^pbf oc v~ x ranges between 22/5 for low-DM pulsars observed at relatively high frequencies (but where 
the scattering is still strong) and 3/2 for high-DM pulsars observed at low frequencies in the super-strong 
scattering regime. This corresponds to the case where the scattered wavefield's coherence length is smaller 
than the inner scale of the wavenumber spectrum. In more detail we have the following cases: Nt <C 1 and 
N v < 1: x = 22/5; N t » 1 and N u < 1: x = 19/5; N t < 1 and N v > 1: x = 11/5; N t > 1 and N„ > 1 
and negligible inner scale: x = 8/5; JV f < 1 and N v <C 1 and where the diffraction scale is smaller than the 
inner scale: x = 3/2. 

Figure [5] shows predicted values of AtiPBF as a function of DM for Galactic coordinates, £ = 45°, b = 0°, and 
particular observational parameters. The NE2001 model for the Galactic electron density (|Cordes &: Laziol 
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2002) was used to calculate Aiigs and Ai/d vs. DM. Different curves correspond to different effective 
velocities that underly the scintillation time scale. The figure shows that the finite-scintle TOA perturbation 
increases rapidly for more distant pulsars and for observations at lower frequencies. The explanation is that 
even though the number of scintles increases rapidly, the pulse broadening time t<z increases even faster. 

4.8. Large Pulse Broadening Times 

For distant pulsars observed at low frequencies and especially those with very narrow intrinsic pulses, the 
pulse broadening time can exceed the pulse width, causing strong distortion of th intrinsic pulse profile 
as well as delaying it. In this case, template fitting must be customized to each pulsar and frequency 
combination. Scintles are too narrow in time and frequency to resolve, so phase retrieva l methods are not 
possible. One approach is to deconvolve the PBF from the measured pulse profile (e.g. iBhat et ail [2003) 
to yield the TOA along with determinations of the relevant PBF shape and the intrinsic pulse shape. This 
requires usage of a family of candidate PBFs for which the deconvolution indicates a best fit that is non- 
negative. The template bank of candidate PBFs needs to include those that are truncated since evidence 
exists, for highly sca ttered pulsars, that scattering screens are spatially truncated transverse to the LOS 
()Lohmer et al.ll2004D . 

Consider a multifrequency timing program where the pulse broadening is larger than the pulse width W 
at lower frequencies and smaller than W at higher frequencies. The purpose of such a program may be 
to conduct precision timing by removing chromatic interstellar effects. The lower frequency measurements 
could be used to identify the form of the PBF, including its characteristic delay, and the frequency scaling 
of the delay. Corrections of the TOAs for the PBF delays and for a contemporaneous dispersion measure 
DM(t) calculated from the TOAs would then follow. This case is discussed explicitly in § [6j 

Another case is where the pulse broadening is large at all frequencies in the program. This situation applies 
to pulsars in the Galactic center because radio scattering is particularly intense so that timing at radio 
frequencies < 12 GHz will not be in the mean-shift regime. TOA shifts from pulse broadening will depend 
in detail on the templates used. In addition, orbital motion that changes the distance to the scattering 
material that encloses the Galactic center will cause secular changes in the pulse broadening. 




t Qms) dm (p c cm ~ 3 ) 



Fig. 5. — Left: Simulated pulse broadening function showing stochastic variations that result from the finite number of 
scintles in the time-frequency plane. This case is for a pulsar scattered into an rms angle of 1.25 mas by a screen 1 kpc 
from both the pulsar and the observer, an observation frequency and bandwidth of 1.4 and 0.1 GHz, respectively, and a 
total observing time of 2 hr. Details of the simulation are available on request. Right: Timing error due to variations in 
the pulse broadening function from the finite number of scintles. Results arc shown for a LOS through the inner Galaxy 
(i = 45°, b = 0°) for two frequencies, 0.3 GHz and 1.4 GHz. with bandwidths of 50 MHz and 200 MHz, respectively. From 
top to bottom in each set, the curves are for effective velocities of 10, 100, 400 and 800 km s . For lower velocities, there are 
fewer scintles averaged over the fixed observing time of T = 300 s so the timing error is larger. The arrows indicate values of 
DM for which the number of scintles sampled N s = 1 . 
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4.9. Angle- of- Arrival and Refractive Scintillation Variations 

Structures in the ISM refract radiation if they are much larger than the Fresnel scale \J \D sa 10 11 cm and 
larger than the so-called refractive or multipath scale -D e ff#d w 1 AU, causing the angle of arrival (AOA) 
to deviate fr om the direct path. Wavefront curvature caused by the same large-scale structures changes the 
flux density ((Rickettl 1 1990) to produce refractive interstellar scintillations (RISS). We give results here only 
for thin screen media, which serve as simple examples and allow scaling laws to be derived. 

Let r be the departure of the AOA from the direct path as viewed by an observer. The induced time delay 
for a fiducial refraction angle of 1 mas is 

At aoa « 7T D <*t% a 1.21 /is # ofr (kpc)0 r 2 (mas) (25) 

where D e s is the effective distance through the ISM. For a thin screen at distance D s from a pulsar at 
distance D from the Earth, D c s = (D — D s ) (D/D s ). (If the refraction angle is defined instead as the 
bending angle produced by the screen, the effective distance is D e g = (D — D S ){D S /D).) The delay scales 
as A 4 if refraction i s produced by a discrete cloud but scale s slightly differently for a medium with a large 
range of scales fe.g. lCordes et alj fl986: Ro mani et alj|1986ft . 

A second effect is associated with the transformation of the TOA to the SSBC ([Foster fc Cordesl Il990) . 
which depends on an assumed direction to the pulsar. Interstellar refraction causes the actual and true 
directions to differ, inducing a perturbation that scales as A 2 , 

A£aoa,ssbc = c n-r® » c r®6 r (t) cos $(t) » 2.4 [is 8 r (mas) cos $(t), (26) 

where Rq/c a 500 s and is the cyclical angle between the LOS n and the Earth-SSBC vector. 

The last effect results from distortion of the scattered image of the source caused by wavefront curvature, 
which modulates the integrated flux by a factor G = G x G y , where G X:V are "gain" variations related to the 
second phase derivatives of the refractive phase in the x,y directions. For simplicity, we assume that the 
scattered image would be circularly symmetric in the absence of wavefront curvature (G x . y — 1), yielding a 

mean delay t^ is °' ) . For G x<y ^ 1 the mean delay becomes 

ipBF,Riss = ^ iso) (G x + G V ), (27) 

as shown in AppendixJEJ Since refractive flux variations are typically tens of percent but can extend to 100% 
or more, this timing modulation can be a large fraction of the mean PBF delay. The variations induced by 
the gains scale as G~ x — 1 oc A 2 . Combined with the A 4 ' 4 scaling of fa, the refraction- induced variations in 
Td are very strongly frequency dependent. Variable p ulse broadening on a time scale of ~ 60 to 200 days 
was reported for the MSP B1937+21 at the 20% level (jCordes et al.lll990t iRamachandran et al.ll2006l) that 
has been attributed to refractive modulation of diffraction, as described here. 

5. A SIGNAL MODEL FOR PULSE ARRIVAL TIMES 

We combine the physical TOA errors highlighted in §§[30 into groups that are convenient for discussing the 
overall timing error budget and for fitting to multi- frequency arrival times. The total TOA perturbation 
has two terms, 

A^toa = Ai w + A£ s i ow , (28) 

one with white-noise statistics and a second that varies slowly enough that it is typically correlated between 
observing epochs. A careful analysis of the errors must also consider the frequency dependences of the 
various terms. The white-noise term is 

Ai w = At j + AtgPBF + Ai s /N . (29) 

The slow term is made up of strongly chromatic contributions, 

A£ s iow = AioM + AtpBF + Ai^pBF : RISS + AtAOA + AiAOA,SSBC + A£dM : i; + Aipjy[- (30) 

The two terms (white and slow) are treated very differently in any arrival time analysis. The white-noise 
terms are unavoidable because they cannot be fitted for and subtracted from arrival times. The slow terms, 
unlike the white terms, can potentially be removed through appropriate multi-frequency fitting because 
there is a finite number of contributing terms that have different frequency scalings. Timing noise and 
gravitational perturbations also produce slow variations but are achromatic. 

We now discuss the white-noise and slowly varying terms in detail. 
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5.1. White Noise Contributions to Timing Errors 

Pulse-phase jitter, the finite-scintle effect and radiometer noise are typically uncorrelated between data 
sets used to estimate individual TOAs. Jitter is weakly chromatic, while the other two contributions are 
strongly chromatic. For radiometer noise, this is due partly to the contribution of Galactic synchrotron noise 
at frequencies less than about 0.5 GHz and to the typical steepness of pulsars' radio spectra. However, when 
interstellar scattering broadens the pulse at low frequencies, the radiometer-noise error is enhanced according 
to Eq. |3] because the pulse- width increases and SNRi decreases. Scintillations are strongly chromatic. 

In Figure [6] we show the three white-noise terms for four different sets of pulse period and DM. The salient 
features of the cases shown in the figure are that 

1. For fixed observing time T, the phase-jitter effect is small for MSPs that have narrow pulse widths 
and for which a TOA is calculated from a large number of pulses. The opposite is true for strong, 
long-period pulsars, which will have jitter dominated TOA errors at frequencies > 0.4 GHz. 

2. The radiometer noise term Aig/N dominates other TOA perturbations for weak pulsars (or 
observations with small telescopes). At low frequencies, the timing error is enhanced by pulse 
broadening from ISS, which also reduces the peak flux density. 

3. For strong pulsars, the TOA error is dominated by radiometer noise at low frequencies because of 
the large Galactic background and by pulse jitter (Aij) at high frequencies. In some cases (e.g. the 
bottom-left panel in Figure [5]) at intermediate frequencies ~ 0.5 — 2 GHz, the finite scintle effect 
(A£pbf) plays a role. 

4. At higher frequencies, the TOA errors of weak, short period pulsars are dominated by radiometer 
noise unless high-gain telescopes (Arecibo, SKA Phase 1, SKA) are used. 

The primary recourse is to reduce the amplitudes of the three contributing terms by increasing the averaging 
time and the bandwidth and by improving the sensitivity of the telescope with a larger collecting area or 
reduced system noise temperature. As a function of frequency At s / N and A^pbf are typically uncorrelated 
(in the latter case for bandpasses that are separated by at least one scintillation bandwidth). The jitter 
term Aij is correlated even between widely spaced frequencies because single pulse structure is typically 
seen over octave bandwidths, though not without some evolution in frequency. A multifrcquency analysis 
may be useful for correcting or reducing pulse phase jitter owing to this partial correlation. 

5.2. Slow, Chromatic Variations 

The slow terms in Eq. [30] are all caused by plasma effects and are therefore strongly chromatic. From left 
to right, the terms scale with frequency approximately as v~ 2 , v~ 4A : v~ x , v~ 2 , v~ 23 / 6 and v~ 3 . We 
emphasize that such scalings are idealized and likely differ between lines of sight and also for particular 
realizations of sampled data. 

In the absence of multipath scattering, the DM term scales as v~ 2 with departures that require very hig h 
electron densities or frequencies near the plasma frequency cx v~ " 4 (|Tanenbaum et alJll968t ICraftl fl970h . 
However, multipath scattering introduces an averaging scale on the DM perturbation that adds the small 
term with i/~ 23 / 6 frequency dependence, as discussed in § 14.41 

The v~ AA scaling of A^pbf assumes a Kolmogorov wavenumber spectrum in the moderate scattering regime. 
In super-strong scattering the exponent becomes —4 because the diffraction scale is much less than the inner 
scale. This regime will apply for high-DM pulsars ob served at low frequencies. Other departures from the 
v~ iA scaling for pulse broadening mentioned earlier (|Lohmer et al.l [20011 12004D indicate that the ISM has 
more complex structure than in idealized models. 

The scaling of the AOA term cx v~ 2 also requires discussion. For discrete structures in the ISM that 
refract the diffracted radiation from a pulsar, this scaling is appropriate. However, the scaling differs for 
density irregularities from a Kolmogorov spectrum because the range of length scales contributing to the 
refraction also varies with frequency. For a Kolmogorov spectrum, refraction effects are small, however, 
so it is appropriate to model the AOA term as we have, given th at there is evidence for the presence of 
non-Kolmogorov refracting structures in the ISM ([Hill et al.ll2005l ). 

The frequency scaling of the RISS modulation of pulse broadening is not a simple power law because 
the dependence on the gains G x<y and on any anisotropy in the diffraction is not simple. However, in 
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the following we assume an approximate scaling law, v XspBF for this term that applies over a restricted 
frequency range. 

The RM term is likely too small to include in any actual fitting model for most lines of sight. It is conceivable 
that a pulsar of interest will be discovered for which profile splitting is important because the RM is large. 
The actual scaling law of the timing perturbation may differ from the generic v~ 3 scaling because any effect 
will depend on how TO As are calculated from the measured pulse profiles in the various Stokes parameters. 

We write terms using values at a fiducial frequency, Vq. In approximate order of decreasing amplitude we 
have, using z = v/vq, 

At s iow(t, v) = AtoMo^ 2 + At A oA,ssBc ^" XAOA + AtpBF ^ XpBF + At A OA ^ 2XAOA 

+At 5 PBF,Riss ^~ X5PBF + At BM ^ z- XDM -" + At nM0 z~ 3 , (31) 

where, for example, for a Kolmogorov medium in strong scattering, Apbf = 4.4, Aaoa = 2, Adm.i/ = 3.8 
and X,5pbf has a number of possible values as described earlier. 

6. MITIGATION OF TIMING ERRORS FROM NON-DISPERSIVE CHROMATIC EFFECTS 

We have identified salient timing errors that are introduced by intervening plasmas. The primary question 
is how to contend with them when precision timing is required, as for the detection of gravitational waves 
with pulsars. We discuss two broad approaches that are based on fitting "raw" arrival times vs. frequency 
v (as opposed to first correcting TOAs for one or more effects with subsequent fitting). 

An augmented version of Eq. [TU]for the TOA includes a non-dispersive, chromatic term, tc, 

U =too +^dmH +tc(v) + t w (v). (32) 

To simplify the notation, we have dropped the 'A' used previously to represent fluctuations. In the following 
we consider multifrequency TOAs at a given epoch. Consequently we include in i-yv only the radiometer 
noise and scintillation effect. The jitter term in this context is considered to be frequency independent (see 
§ I3.3[) . When using the variance of the white noise in weighted least-squares fitting across frequency, we 
exclude the jitter term for this reason. However, the jitter term does contribute to epoch-to-epoch timing 
errors, as previously discussed. 

6.1. Ignore Refraction and Scattering Effects 

One approach is to simply ignore all effects other than the dispersion delay. Most arrival-time analyses 
to date include the dispersive term A^dm as the only chromatic perturbation. Such programs use dual or 
multiple-frequency TOAs to estimate DM at each epoch or use ancillary measurements of DM to correct 
TOAs to infinite frequency a t the SSBC. This procedure allows scattering to contaminate the estimated 
TOA ([Foster fc Co rdes 199Cj). To be successful, removal of dispersion delays only requires additional errors 
to be smaller than some target precision. As in Figure HI for the case of pulse broadening, restrictions can 
then be made on the combination of observation frequencies and DMs of a pulsar sample that can satisfy 
the error budget. 

Following the same approach discussed in § 14.11 and in Appendix [D] we calculate the error on the corrected 
TOA, too, when the tc term is present in the data but only the £dm term is included in the fit. For example 
we consider the simple case where that tc includes only one term instead of the many terms of Eq. [31] and 
that it scales as tc{v) — acv~ x , (where we drop any subscript on the exponent, X). The two parameter 
fit for too and DM yields a systematic error for too ( m addition to the random, standard error) that can be 
written as 

Stoo = -a c R toa , (33) 

where Rt^ is a coefficient defined in Eq. ID3I For X — 4, a 2:1 range of frequencies, and frequency channels 
ranging in number from 2 to 1024, we find a variation i? too « 0.3 ± 0.05. Thus a fraction 1 — i? too ~ 70% of 
the non-dispersive delay ac is absorbed in the fit for DM, yielding an estimate for too that is too early by 
~ (0.3 ± 0.05)ac. For a 4:1 frequency range, R tx ~ —0.1, so 90% of ac is absorbed. However, for a 1.5:1 
range, only about 50% is absorbed in the fit. Referring to Figure|3]for a 2:1 frequency range, the maximum 
timing error at from this effect alone requires that ac < o't/Rt^ or, since ac = tc{v\)v x , 

tciyx) < vI x R7^t- (34) 

For vi — 2 GHz, <r t = 0.1 fis and using i? too = 0.3, we need ic(2 GHz) < 19 ns. 

A more detailed analysis is shown in Figure [7] with calculated results for an MSP with P = 1.5 ms with low 
DM (10 pc cm" 3 ) with Galactic coordinates £, b — 45°, 0° and a pulse width of 50 /is. The period-averaged 
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flux density is 1 mJy at 1.4 GHz and scales with frequency as v~ 2 . The NE2001 model is used to estimate 
the distance and pulse broadening from scattering. In the left hand panel, the rms timing error is plotted 
against bandwidth B for the case of an MSP observed with a 100-m class telescope. The right-hand panel 
shows the corresponding errors for DM. The example shown in the figure applies to observations using 100- 
m class telescopes. The highest frequency in the example is 2 GHz so the measurements and fit therefore 
cover the frequency range [2 — B, 2] GHz for each bandwidth plotted. The left-hand panel shows standard 
errors for fitting functions with one to three parameters (too, DM an d ac) and systematic errors for the 
three cases. Other details are described in the caption. The two-parameter fit is relevant to this section 
where pulse broadening is ignored. Salient features of the figure include: 

1. For the k = 2 fit for too and DM the total error minimizes at B 1.4 GHz as pulse broadening 
becomes important at the lower frequencies in the band. 

2. When pulse broadening is ignored in any fitting, the best approach is to use a two-parameter fit if 
DM is not known to sufficent precision to warrant a single-parameter fit. 

3. For small bandwidths, standard errors (dashed lines labeled csE.k, k = 1,2,3) become increasingly 
larger in going from one to three-parameter fits because the weak dispersive and scattering terms 
are not well modeled. The timing error scales as i?™ 1 / 2 , _B~ 3 / 2 and B~ 5 ^ 2 , respectively, when the 
TOA error is dominated by radiometer noise, as it is for small bandwidths. For B -C V\ the ratio 
cse,2/cse,i ~ y/Zvi/B and is a factor of 69 for the smallest B plotted in Figure[7] The penalties for 
fitting for DM and to can be avoided if the dispersion measure and pulse broadening are known to 
sufficient accuracy from ancillary measurements, but both of these quantities vary with epoch and 
must be monitored. 

4. Foreshadowing results in the next section, the three parameter fit improves the TOA error but 
only if the scaling exponent A is known. The upturn in or,3 is caused by a 10% uncertainty in A, 
SX = 0.44, illustrating the requirements on prior knowledge of A. Nonetheless the improvement is 
by a factor of about 10 and can be greater if A is known to better precision. 

5. The best achievable timing precision corresponds to (Jt,s and is about 0.1 [is for the 600-s 
observation shown and with a 1.8 GHz bandwidth. 

6. Brighter pulsars will shift the curves downward over most, but not all of their plotted range. 

We emphasize that the results shown in Figure [7] and in the next two figures do not include all refraction 
and diffraction effects, so they must be taken as illustrative and optimistic for considerations of ISM effects 
alone. 

6.2. Aggressive Removal of Refraction and Scattering Effects 

If non-dispersive effects arc tackled more aggressively, a larger number of pulsars can be included in a timing 
program because a wider range of frequency-DM combinations can satisfy a given error budget. 

Consider the case where, along with too and 4dm j we include tc in the fitting function as a power-law 
tc = a c v~ X ■ 

This approach does not require knowledge of C\ or measurement of the scintillation bandwidth Ai/d as 
discussed in § |4.6.1l but it does require assumption of a value for the exponent, A. For now, we consider A 
to be known and we consider tc to describe only pulse broadening. 

Figure [8] shows results similar to those in Figure [7] for observations with v\ = 2 GHz of a 1.5 ms pulsar at 
DM = 100 pc cm -3 using a 100-m class telescope (left panel) and a 300-m class telescope. The latter case 
corresponds approximately to the Arecibo telescope, the Chinese FAST telescope, and the proposed Phase 
1 of the SKA. 

The standard errors shown in the figure are larger for the 100-m class telescope in accord with its smaller 
sensitivity by about a factor of four. The total error from the two-parameter fit minimizes for B « 0.5 GHz 
and B ~ 0.4 GHz for the smaller and larger telescopes, respectively, the upturn being caused if there is 
no correction for pulse broadening. A similar upturn occurs for the three parameter fit if there is an error 
in the scaling-law exponent 8X — 0.44 (10%). If there is no error, the best rms error is for a bandwidth 
B = 1.6 GHz and is approximately 2.5 times better for the larger telescope. 

The conclusions that can be made from this particular case are that aggressive treatment of pulse broadening 
can drastically improve timing precision of moderate-DM pulsars and thus allow more pulsars to be included 
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in a pulsar timing array than otherwise would be the case. Without addressing pulse broadening, there is 
a maximum usable bandwidth that is less than 1 GHz and the timing precision ~ 1 (is is not competitive 
for PTAs unless more observations are made during a year. 

Figure IH1 compares observations of a high DM pulsar (500 pc cm -3 ) at high frequency (5 GHz) for obser- 
vations with an Arecibo-class instrument (left panel) and an SKA-class instrument (right panel). For the 
SKA we have used a nominal sensitivity A c g/T sys = 10 4 m 2 K _1 . Note that the spin period in this case 
is 3 ms. The results indicate that aggressive fitting for pulse broadening can yield timing precision that is 
competitive for PTA applications even for a high-DM object. 

The examples shown in Figures [7] - [9] that include only pulse broadening can be used as a guide for how to 
deal with the other chromatic terms in Eq. 1311 In total there are eight unknowns to be solved for at each 
epoch when we include the achromatic term t^. However, the RM term can almost always be ignored. The 
term A£aoa,ssbc is or almost is degenerate with the DM term in its frequency dependence oc v~ 2 and three 
of the other chromatic terms (AipBFj A^aoa and AiDM,i/) are degenerate with scalings close to v~ A . This 
suggests a minimum of four unknowns that require four measurements at each epoch over a frequency range 
that is adequate to solve for the unknowns. Additional terms can be fitted for to potentially improve the 
net timing precision. However, examples shown in the figures indicate that the "cost" of extra parameters is 
high and that the amplitudes of additional effects need to be large enough to surpass the cost in precision. 

6.3. Comparison with Well-timed MSPs 

We compare our predictions with actual measurements on four prominent MSPs to further illu strate timing 
budg et considerations. All four objects show TOA precision of 200 to 300 ns in 5 min averages (|Hotan et al.l 
2006) where, to be explicit, we mean the precision of the TOA estimate and not necessarily the rms residual 
over a time span of years. Three objects (J0437— 4715, J1713+0747, and J1909— 3744) show rms timing 
residuals that, when integrated over ~ 1 hr, are below 100 ns. The fourth, B1937+21 (J1939+2134), is 
bright but shows timing noise consistent with that seen in canonical puls ars, which is caused by torque 
variations. We use flux densities, pulse widths and other parameters from iHotan et al.l (|2006l ) and other 
references, as noted, an d use parame t ers fo r the Parkes telescope and the 48 MHz bandwidth and 300-s 
integration time used by IHotan et al.l (|2006| ). 

J0437— 4715 is the brightest of the four pulsars, with a period-averaged flux density S1400 = 137 Jy. Eq. [3] 
combined with the pulse width W = 85 (is and the system noise of the Parkes telescope (S sys = 29 Jy) 
predicts At s / N = 13 ns and Atj = 105 ns. This object, not surprisingly, is very heavily jitter dominated. 
The low-DM (2.6 pc cm~ 3 ) implies that interstellar effects will be much smaller than either of these leading 
effects. Increases in sensitivity with larger bandwidths or larger telescopes will not improve the TOA 
precision, but longer integration times and a larger numbe of TOAs will. 

J1713+0747 is less bright and has a larger DM (16 pc cm 3 ) than J0437— 4715. For the Parkes telescope, 
the radiometer noise TOA error is about 50% larger than the jitter error (100 ns vs. 65 ns), but for larger 
telescopes jitter will dominate. Diffractive interstellar scintillations will strongly modulate the signal level 
at most frequencies, making even the Parkes measurements jitter dominated some of the time. 

J1909— 3744 is one of the best timing pulsars. The low DM (10.4 pc cm~ 3 ) is similar to that shown in 
Figure [7] and implies that interstellar effects are small at frequencies above 1 GHz. The measured rms 
residual is a t — 0.3 (is. We estimate Ais/N ~ 0.41 (is, consistent with the curve labeled 'S/N' in Figure[7] 
(left-hand panel) for the bandwidth used and after taking into account the small differences in values for 
period, width, flux density and system noise as well as the radio frequen cy and spectra l index . iBailesI (|2007|) 
notes that the pulsar flux density varies from DISS by a factor of 30 and l Jacobv et al.l (|2005l ) exclude TOAs 
with errors in excess of 1 (is. The difference between the predicted and actual rms residual is attributable 
to modulations by DISS and by RISS. The narrow pulse for this object (42 (is) implies that pulse jitter 
does not make a significant contribution to the observed measurements (Aij « 27 ns from Eq. [5]) unless 
the scintillated flux becomes very large. 

The fourth object, B1937+21, has the shortest period (~ 1.6 ms) and narrowest pulse (30 (is) of the four 
objects we consider. It also has the largest DM (71 pc cm -3 ) and is fairly bright. For the Parkes telescope, 
the TOAs are noise dominated with estimates Atg/N — 56 ns and Aij = 19 ns but will be jitter dominated 
for most larger telescopes. As noted earlier (§ 13. 2\i . jitter in this object is ambiguous, displayed in giant 
pulses but not in typical pulses, so our estimate of Atj is uncertain and there is the possibility, not yet 
demonstrated, that timing precision can be improved by excluding giant pulses from timing estimates or by 
otherwise correcting for them. 
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7. SUMMARY AND CONCLUSIONS 

Table [2] summarizes how timing precision may be limited by different causes and what the course of mitiga- 
tion might be. We include timing noise intrinsic to the pulsar because many pulsars w ill in fact haye tim ing 
residuals that are dominated by timing noise unless it can somehow be removed (e.g. lLvne et al."ll2010l ). 

For many objects, timing precision will be limited by a combination of radiometer noise and pulse-to-pulse 
phase jitter intrinsic to the pulsar. We have summarized the levels of jitter in different pulsars as presented 
in the literature and we have presented specific results on the bright MSP, J1713+0747. We have also shown 
the relationship between jitter and pulse-shape stability. The theoretical best timing precision is achieved 
when it is limited by additive white noise alone and it is straightforward to calculate the rms timing error 
Ai s / N using a pulsar's template pulse shape f Appendix |A"]) . Comparison of actual TOA errors against Atg/ N 
shows that for all pulsars except the MSP B1937+21 (|Jenet & Gilll2004f ). the rms phase jitter appears to 
be a significant fraction of the pulse width and thus will limit timing when the signal-to-noise ratio is high. 

In the noise dominated regime, TOA precision may be improved by actions that favorably alter any factor in 
Eq. lAll increasing the sensitivity (i.e. reducing S sys ) and the total bandwidth to increase the signal-to- noise 
ratio; improving dispersion removal to minimize the effective pulse width W; and increasing the number of 
pulses summed. In the jitter dominated regime, however, TOA precision is improved only by summing more 
pulses. Thus in either regime, an increase in telescope time per object is needed to decrease the timing 
errors from radiometer noise or jitter. As more sensitive instruments are built, more pulsars will have 
jitter-dominated timing, as we have shown in § 13.31 With the Square Kilometer Array, for example, most 
of the known canonical pulsars and the majority of MSPs will have jitter-dominated TOAs. Consequently, 
the SKA can be used with sub-arrays in order to increase the amount of observing time per object with the 
sensitivity of a sub-array set so that the TOAs are jitter dominated by only a small factor. 

Timing variations from interstellar dispersion and scattering need to be characterized for each LOS owing 
to the richness of structure in the ISM. Low-DM pulsars (e.g. < 50 pc cm -3 ) are clearly better for precision 
timing at frequencies ~ 1 GHz but the maximum DM is smaller for lower frequencies. For pulsar-timing 
array applications that use an ensemble of MSPs, the best ap proach is to find spin-stab le MSPs with small 
dispersion measures. However, intrinsic spin noise in MSPs ((Shannon fc Cordes Il2010h may require usage 
of highcr-DM MSPs and mitigation of interstellar effects. We have shown that multiple frequency fitting of 
TOAs obtained over a large total range of frequencies can extend the allowed sample of pulsars. Objects 
in relativistic binaries are rare enough that future discoveries are likely to include pulsars with large DMs. 
The most extreme case is for pulsars orbiting the Milky Way's central blac k hole, Sgr A*, whic h will show 
extreme levels of scattering with corresponding large timing variations fe.g. lDeneva et alJl2009t ). 

The dominant interstellar timing error is from variations in DM followed by diffractive pulse broadening 
and its variations. Both of these effects are correctable to vary degree. Additional interstellar terms imply 
that timing measurements at a minimum of four narrowband frequencies over approximately a 2:1 frequency 
range. Better yet — and now feasible — are continuous measurements over octave-like frequency ranges. 
For small dispersion measures, the pulse broadening time is too small to be determined from TOAs alone 
even though it is large enough to influence the arrival times. We have outlined a procedure that fits nominal 
TOAs at multiple frequencies to remove chromatic effects. However, other approaches, also mentioned here, 
could alternatively remove some of the chromatic effects at individual frequencies before doing any multiple 
frequency fits. Further exploration is needed on using dy namic spectra and their analy sis using either 
correlation functions or a secondary spectrum analysis (e.g. iHemberger fc Stinebring||2008f) to provide the 



Table 2 

Regimes for Pulsar Timing Limitations 



Limiting Effect 


Spectral Signature 
in TOAs 


Regime Specifications 


Mitigation 


Timing noise 


red 


spin, torque noise 


more pulsars in timing array, 
remove torque variations 


Low S/N 


white 


weak pulsar, small telescope 


larger T, B\ smaller S ByB 


Phase jitter 


white 


large S/N 


larger T 


Diffractive scintillations 


white 


moderate DM 


large T,B,v 


(Finite scintle effects) 






S/N + pulse broadening 


white 


large scattering 


larger v or reject pulsar 


DM, AOA variations 


red 


moderate DM objects 


multi-frequency fitting 


PBF + variations 


red 


slow changes in t& 


multi-frequency fitting 
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necessary timing corrections. Other approaches that use phase information in high resolution sampling (e.g. 
P. Demorest, W. van Straten & M. Walker, private communication) also need to be explored. 

We thank P. Demorest, B. Rickett, D. Stinebring, and M. Walker for useful discussions or correspondence 
and the Nanograv collaboration and the Lorentz Center for providing useful venues for discussions of timing 
precision. JMC had very fruitful and enjoyable conversations about timing precision with Don Backer 
over many years that have contributed to this paper. This work was supported by the National Science 
Foundation, which supports the Arecibo Observatory under a cooperative agreement with Cornell University. 
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A. TIME-OF- ARRIVAL ERRORS FROM PULSE PROFILE STOCHASTICITY 

Here we give expressions for the error in time of arrival (TOA) due to additive radiometer noise and to 
pulse-phase jitter that is intrinsic to the pulsar. 

Radiometer noise is additive in the limit of large time-bandwidth product and causes a minimum error 
when matched filtering is used to estimate the TOA. Let U(i) be the pulse waveform as a function of time 
and normalized to unity amplitude. Also define the autocovariance function of the noise to be p(r) = 
(n(t)n(t + T))ja\ oise such that p(0) = 1 and its characteristic wi dth is A = [ dr p(r). Fo r a signal-to-noise 
ratio SNR (pulse maximum/rms noise), the rms timing error is (jDowns fc Reichlevlll98l 

a , _ [fJdtdt'p(t-t')U'(t)U'(t')] 1/2 _ W cS (A\ 1/2 = S svs [W cB \ 1/2 
S/N SNR/dt[U'(t)] 2 SNRxVNVWW S poak V2B\Nj ' [ ' 

where U' = dlf/dt. The second equality holds in the limit where noise decorrelates on a time scale A 
much less than the pulse width, i.e. A <C W for which p(t) = AS(t) and where the effective pulse width is 
W e B = 1// dt[U'(t)} 2 ; the result is written in terms of the SNR for a single pulse and the number of pulses 
summed, SNRi = SNR/vN- The third form shows that the error does not depend on A, where we have 
used S sys as the system equivalent flux density, S^eak as the average peak flux density of a pulse, and B is 
the total bandwidth. The factor of two accounts for summing of two polarization channels. 

For measured pulse profiles Uj — U(jAt) sampled discretely at intervals At <C W a s, the effective width can 
be calculated by differencing, 

W cS = — — -. (A2) 

£([/ i+1 -^-) 2 

3 

The denominator is proportional to the structure function of the pulse profile. If the sampled pulse has 
finite signal-to-noise ratio, the denominator must be corrected for the additive noise by subtracting 2n on <7 2 , 
where a n is the rms of the additive noise and n on is the number of terms in the sum. For a Gaussian pulse 
with width W (FWHM), 

A<S/N = (2^1n2) 1 /4 S NR 1 VN \ W ) ' (A3) 

We describe jitter by letting a(cf>) be the functional form of a single pulse and f<t,(4>) the probability density 
function (PDF) for the phase pulse-to-pulse phase jitter. Assuming all pulses have the same shape, the 
ensemble-average pulse shape is 

U(<j>)<x J dtf U{<t>>{4> ~4>')- (A4) 

With this model and the phenomenologically consistent assumption that an infinite pulse average converges 
to U(4>), the TOA error is (using the pulse period P) 



At, = n; 1/2 {\ + m 2 ) 1 / 2 p{<t> 2 ) 1/2 - n; 1/2 {\ + m 2 ) 1 ^ P 



d^4> 2 U(^ 



1/2 



(A5) 



where Ni < N is the number of statistically independent pulses that have been summed and mi is the 
intensity modulation index (ratio of rms intensity to mean intensity), typically of order unity. Intensity 
fluctuations without phase jitter yield no timing error but they enhance timing errors when there is phase 
jitter. Letting (0 2 ) 1 / 2 = £,ufjWi/P we define a dimensionless jitter parameter fj equal to the width 
(FWHM) of the phase-jitter PDF in units of the resulting intrinsic pulse width Wi, where £u ~ 1 is a 
shape-dependent factor. For Gaussian shaped pulses with a Gaussian phase jitter distribution, we obtain 

Atl ~ 2(27^2)1/2 • ( A6 ) 

This idealized result will also apply approximately to pulsars with multicomponent profiles. The amplitudes 
and phases of individual subpulses in different components are only partially correlated, so a profile with 
M components would show a jitter-induced error that is up to a factor M" 1 / 2 smaller than that for a single 
component. 
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B. SCATTERING AND REFRACTION FROM THIN SCREENS 



Electron-density fluctuations on large scales refract, focus and defocus radiation from a pulsar producing 
variations in flux density and arrival time. "Large" is denned relative to the scale ~ D9 c i that is larger than 
the Fresnel scale ~ V XD by the same factor that the Fresnel scale is larger than the diffraction scale ~ X9d- 
(e.g. Rickett f990). While an acceptable model for electron-density variations is one with a continuum of 
length scales, it is useful to separate large and small scales, as we do here. The Kirchoff diffraction integral 
for the scalar wavefield has an integrand exp [i$(x, x')], where x and x' are both two dimensional vectors 
transverse to the LOS. The total phase $ = 4> g + <p r + <pd comprises a geometric term 

j 2 



; (x,x)=(- 



D~ x x 



(D - A 



rV-x'i 2 



(Bl) 



a refraction term, <j> r , from large-scale variations and the diffraction term, fa, from small-scale variations. D 
is the Earth-pulsar distance and D s is the pulsar-screen distance. In the absence of the refraction term, the 
scattered wavefield of a point source is described by an ensemble-average image Io(6). We include refraction 
as a quadratic surface centered on a stationary-phase point x' on the screen corresponding to an angle 8, 

</> r (x') = fa(x') + b ■ (x' - x') + (x' - x') • C • (x - x'). (B2) 

The stationary phase point is defined with respect to fa+fa and generally there could be multiple stationary 
phase points. By consolidating linear and quadratic terms, the refraction-distorted image can be written 

as 



where the gains are given by 



I(6) = I (G- 1 (9 X -9 X ),G- 



/)). 



(B3) 



G x . y = [1 + k- L (D s /D)(D - D s )C Xt y\ . (B4) 

We assume, for simplicity, that ellipses of constant phase from the last term in Eq. IB2I align with the 
transverse coordinate axes (x,y). In this case, C, is diagonal with elements C X:V . A gain G XiV > 1 implies 
that the image is wider in the x, y directions compared to the unity gain case. While the amplitude of the 
image is unaltered, the integral over the image implies a flux variation by a factor G x G y . A more general 
expression that we use in Paper III is 

/(e) = /oCu t Ts- 1 u(e-e)), (B5) 

where U is a 2 x 2 matrix that diagonalizes C, the dagger denotes transpose, and 72 is a diagonal 2x2 
matrix with elements G" 1 and G y . 

For each point on the screen the time delay relative to the direct ray in vacuum is 



1 

t — tr>u + 7T 

2c 



D{D-Ds) 



9 +G- 



G: 



(B6) 



The first term £dm = (27r^) _1 (/)o(5f') is the extra dispersive delay from the position of the image centroid. 
The second term is the geometric delay from the centroid of the refracted image. The third term results 
from the combination of extra geometrical delay and modified dispersive delay from the quadratic part 
of the refractive phase. A gain G x > 1 requires a converging wavefront that is produced by a deficit of 
electrons compared to the local mean in the screen. Thus while G x > 1 increases the size of the image in the 
x direction and increases geometric delays oc G x , these are compensated in part by the reduced dispersion 
delay. 

We average the last term over the distorted image, which we treat as a probability density function with 
angular variances (G x a x0 ) 2 and (G y a yo ) 2 , yielding 



1 

t — £dM + 7T 

2c 



D(D — D s ) ~ 

When the undistorted image is symmetric we have a x0 
is the the sum of the gains, G x + G y . 



-2 



9 +G X&XQ 



G y a yo 



(B7) 



a yo and the time delay then involves a factor that 



C. FREQUENCY DEPENDENT DISPERSION MEASURE FROM PHASE-SCREEN AVERAGING 



Scattered pulsar radiation reaches the Earth from a cone of radius 9^, the diffaction angle. For simplicity 
here we assume scattering occurs in a screen a distance D 1 = D — D s from the Earth. The screen phase 
</>(x) is averaged by the area of the scattering cone ~ t 2 s where l s — D'9d- DM is related to the screen phase 
by DM = ^0/27raDM where cjdm is defined in Eq. [9j The screen-averaged DM is 



DM(z/) 



dxW(x-x ,z/)DM(x), 



(CI) 
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where the averaging function W(x) has unit area and has a width that is strongly frequency dependent. 
Letting Xo be the point on the screen where the direct ray path emerges with an associated dispersion 
measure DMo = DM(xq), we derive the variance of the difference <5DM(V) = DM(V) — DMo. Assuming that 
the screen phase has stationary statistics, that scattering is circularly symmetric, and assuming that the 
weighting function has an autocorrelation function i?w(x, x', v) that depends only on the difference x — x', 
we can write the DM variations in terms of the phase structure function D^(5x) 

2 



([SBM(u)} 2 ) 



1 
2^ 



«DM 



dyyD^y) 



1 



W(y,v) - -R w {y,v) 



(C2) 



If we use the Kolmogorov form D^y) oc y 5 / 3 (see main text and lCordes fe Laziol (|1991h ; [Cordes fe Rickettj 
(1998)) and if we approximate the weighting function as £j 2 that drops to zero at a radius £ s , we obtain 



([5DM(v)f) 



3/, 



5/3 



22tt 2 



"DM 



SM£ 5 J 3 



(C3) 



Using i s = D'Bd and the relation for a Kolmogorov medium in strong scattering (|Cordes fe Laziol I1991L 
2002) 



6 d = 0.071 arc sec SM 3/5 ^~ 11/5 , 

we obtain 

([<SDMM] 2 ) 1/2 
The associated rms time delay is 



10" 1 ' 05 pc cm" 3 SMD 



/5/5 -11/6 



AtsDM = 4.15 x 10" 3 s v- 2 ([<5DM(i/)] 2 ) 1/2 w 0.118 fxsD' 



/5/6 -23/16 



SM 

IF 5 * 



kpc m 



-20/3 



(C4) 
(C5) 

(C6) 



D. LEAST-SQUARES FITTING OF MULTIFREQUENCY DATA 

Consider a set of multi-frequency TOA measurements obtained in an observing session at a given epoch, 
{tk, Vk, k = 1, n^}. We model the TOAs as t(y) — t o+^DM(^)+ic( ;y ) + e ( i/ ) where is the TOA at infinite 
frequency, £dm is the delay from plasma dispersion, tc is an additional, non-dispersive chromatic pertur- 
bation (e.g. from scattering), and e is measurement error that we consider to be white noise (statistically 
independent). Actually, tc is the sum of several effects as outlined in the main text, but for convenience 
here we use only one term. In vector form the TOAs are D = X0 + e, where D = col{tk} is a column vector 
of arrival times, is a vector of parameters, and X comprises the independent variables. With frequency- 
dependent errors described by a covariance matrix C, the solution is — (XtC _1 X) X^C X D. The 
covariance matrix for parameter errors is P = (X^C _1 X) 1 if the model is correct (only additive errors 
to the TOAs). For white-noise errors in TOAs at the different frequencies, C = diagonal {c^} ■ ^ n the 
following we will use the quantities 

r p = ^o-fe 2 ^ p ln^, s p = ^o^ 2 v- p . (Dl) 

fe k 

If we fit only for and there are no chromatic effects (the situation at very high radio frequencies and 
in high-energy observations) the variance is af = 1/sq. When the errors are identical at all frequencies, 

of^ = <j 2 1 fly. 

Fit for dispersion: Consider a two-parameter fit that includes the dispersion term along with and 
that there are no other chromatic effects in the data. The parameter vector is 9 = col (too, «dmDM), the 
variables matrix is X = matrix(l i/r ), i = 1, . . . , n„, and the standard errors are of = S4/D2 for and 
"dm = s o/flDM^2 for DM, where T>2 = S0S4 — sf, is the determinant of X^C _1 X. 

An alternative approach uses an ancillary value of DM (e.g. interpolated from measurements at other 
epochs) that will have an error <5DM. The TOA error is then a 2 ^ = Sq 1 + (aDM<5DMs2/so) , which is better 

than the two-parameter fit if <5DM < aj^ (so/2?2) Evaluating in the continuous limit with errors that 
are independent of frequency and letting = v\ — B where B is the total bandwidth and <jb = o 7 1 ^fnZ^ 
we require 

5DM < aB ^ (1 - B/^) 3 / 2 = 10- 3 - 38 pc cm" 3 ^(1 - B/^) 3 / 2 a B M (D2) 
for vi and B in GHz. 
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Now consider the case where the data contain the chromatic term tc(y) but the fitting function does not. 
If tc{v) is deterministic^ the additional, systematic errors are 

St x = V 2 X ^2 a k 2 ( S4 ~ s 2vl 2 )tc{v k ) = ~a c R toa = acD^ 1 (siSx - s 2 s x+2 ) (D3) 

k 

SDM = V^ 1 ^2 a k 2 ( s o l/ k 2 ~ s 2)tcOk) = a c (a DM X>2) _1 (s s x +2 - s 2 s x ) (D4) 
k 

where the last equality in each equation holds for a power law scaling, tc{v) — a c v~ x . Because it absorbs 
some of the scattering delay, tc, SDM is always positive, causing TOAs to be overcorrected and producing 
a negative bias in too. 

A simple case is where narrowband measurements are made at two frequencies v\$, with r — v\/v2- The 
error from the chromatic term is St^ = (f 2 tci — tc 2 ) I if 1 ~~ l) ■ F° r an octave bandwidth (r = 2), the TOA 
error is approximately —1/3 of the chromatic term at the lower frequency when v\ 3> V2 so that tc 2 3* tci ■ 
Not surprisingly, observations over an octave at low frequencies where tc is larger will be substantially 
worse than at high frequencies. For large r we have St^ — r z tc i: showing that the chromatic term at 
the higher frequency is amplified by a large factor when X = 4. 

For n v — 1024, v\ = 2 GHz and v\ = 1 GHz, we get St^/dc = —0.29 and 5DMaDivi/ac = 1.15. Suppose 
that ac = 1 (J-s if all frequencies are expressed in GHz, equivalent to fc(l GHz) = 1 fis. In this case 
^oo = —0.29 fis and SDM — (1.15 /zs/4.15 ms)pc cm~ 3 = 10~ 3 ' 62 pc cm~ 3 . It may be more convenient 
to express the systematic errors in terms of tc at either the highest or lowest frequency, in which case the 
relations ac — tc v \^ = tc v 2 can be used. 

Various trends can be identified by considering other cases. For fixed total bandwidth B, the systematic 
TOA error decreases with increasing upper frequency V\. For fixed upper frequency, the systematic error 
increases as B gets larger (i.e. as the lower frequency gets lower). Counteracting this trend is the decrease 
in standard error with larger bandwidth, implying that there is an optimal bandwidth that depends on the 
pulsar spectrum and system noise vs. frequency along with the scattering strength. 

Fit for dispersion and scattering: Now we include scattering as a single term in the fitting function, 
tc{y) — clcV ~ i and assume initially that the exponent X is known. Then — col (too, odmDM, ac) and 
the variables matrix is X = matrix(l v^ 2 v^ x \ i = 1, . . . , n v . 

When X is assumed to be Xq, the error SX — Xq — X, causes a systematic error from an extra term 
obtained from the expansion, tc(v, Xq) w ac^ x — oq5Xv~ x In v. The error on the arrival time to linear 
order in SX is 

St x = ~a c SXG too (D5) 

where 

G *=o = T^z 1 [ r x (saS2X - 4+2) r x+2 {sxsx+2 - s 2 S2x) + r 2 x (s 2 s x +2 - S4,s x )] (D6) 

and 

X> 3 = dot X t C _1 X = (s Si - s 2 2 )s2x + 2s 2 s x s x +2 - (s4,s 2 x + s Q s 2 x+2 ). (D7) 

Evaluting for the same case as above with n v — 1024, v\ = 2 GHz and Vi = 1 GHz, we obtain Gt x — —0.069 
implying for SX > that, opposite to before, the TOA is underrcorrected because DM is biased low while 
the scattering parameter ac is biased high. When SX < the opposite is true. Either way, the amplitude 
of the timing error is small, Stoo = 0.0Q9acSX . 

Finally, a four parameter fit can include SX in an iterative approach that uses an initial value for X 
and interates. The parameter vector is = col (too, odmDM, ac, ac <5X) and the variables matrix is X = 
matrix(l v~ 2 v^ x — ac v 7 In ^i), i = 1, . . . , n v . In this case, initial choices for oq and Xq are updated 
by fitted values for ac and SX. 



3 

Stochastic tc will yield an error in ioo that involves the correlation function {tc{v)tc(v'))> but we will not elaborate on this 
case in this paper. 
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l.b=30., 0. G= 1. 70000005 K Jy' T jn[ =600. s l,b=30., 0. G= 1.70000005 K Jy- 1 T lnl =600. s 

P-0.0015 DM=10 pc cm" 3 S 110D = 1 000 mjy spectral index = 15 P=0 0300 DM=10 pc cm- 3 S 1400 = 10 000 mjy spectral index = 15 




Frequency (GHz) Frequency (GHz) 

Fig. 6. — Rms timing errors that contribute to the white-noise term in Eg. 1291 are shown plotted against frequency. The 
top two curves are calculated for a 100-m class (GBT-like) telescope with gain G = 1.7 K Jy — 1 while the bottom curves 
are for an Arecibo-like telescope with large gain, G = 10 K Jy — . Flux densities quoted below are for 1.4 GHz and scale 
with frequency as u~ 15 . The bandwidth ratio B/v is 0.05 for frequencies below 0.3 GHz, increases linearly to 0.2 from 
0.3 to 1 GHz, and stays at 0.2 for higher frequencies. Solid (red) curves: the total white-noise contribution to the TOA 
error. Dot-dash curves: Atg/ N calculated using the intrinsic pulse width; these are underestimates of the true TOA error 
from radiometer noise. R adiometer noise includes a constant temperature added to the 408 MHz background temperature 
from JHaslam et al.|[l982l) . The intrinsic pulse shape is Gaussian in form with a duty cycle of 0.03 at 1.4 GHz that scales 
as u~ . The duty cycle also scales with spin period as P~ 0,4 for periods longer than 2 ms but are allowed to be small 
(0.02) for shorter periods to model those MSPs with very narrow pulses. The pulse width used to calculate the contribution 
from radiometer noise (Atg/ N ) includes the effects of interstellar pulse broadening calculated using the NE2001 electron- 
density model. Dotted curves: the radiometer-noise timing error (Atg/ N ) calculated using the apparent pulse width that 
includes interstellar scattering broadening; Short-dashed curves: The contribution from pulse-phase-jitter (Atj) using a jitter 
parameter fj = 1/3. Long-dashed curves: Timing errors from the finite number of scintles (A^pbp). The scintillation time 
used to calculate At^PBF is based on a fiducial transverse velocity of 100 km s . (Top left:) Millisecond pulsar (P = 1.5 ms) 
that is nearby with DM=10 pc cm -3 ; the TOA is noise limited because the flux density (1 mjy) and the telescope gain 
(1.7 K Jy — x ) are low. (Top right:) 30-ms pulsar that is nearby but bright enough that its timing is jitter dominated for most 
of the frequency range. (Bottom left:) Millisecond pulsar with moderate DM = 100 pc cm -3 and high enough flux density 
(10 mjy) observed with high gain (10 K Jy — x ) so that the timing is dominated by the DISS finite-scintle effect. (Bottom 
right:) 30-ms pulsar with large DM (300 pc cm~ 3 and large flux density (10 mjy) observed with high gain such that the 
timing errors are jitter dominated at high frequencies but are noise limited at low frequencies owing to pulse broadening from 
interstellar scattering. 
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Fig. 7. — Left: Standard and systematic errors after fitting multifrequency TOAs across the designated bandwidth with 
a maximum frequency of 1.5 GHz. These errors apply to cases where the only chromatic terms included are the dispersive 
term *dm and an additional term tc oc u~ x with X = 4.4 that approximates pulse broadening or one of the refraction terms 
discussed in the text. The curves apply to a 600 s observation using a 100-m class telescope of a pulsar with P = 1.5 ms and 
DM = 10 pc cm' 3 in the direction I = 45°, b = 0°. The pulse width is 50 [is. Dashed lines labeled <rsE,ki fc = 1,2, 3 are the 
standard errors for one, two and three-parameter fits as discussed in the text. Solid lines labeled ^ are total errors that 
quadratically sum standard and systematic errors. For k = 1, the systematic error includes an error in dispersion measure, 
<5DM = 10 — 3 5 pc cm -3 along with pulse broadening calculated using NE2001. For k = 2, the systematic error is solely 
from pulse broadening because DM is fitted for. For k = 3, an error SX = 0.44 (10 %) is assumed in the exponent of the 
pulse- broadening scaling law. The curve labeled "S/N" is the standard error resulting solely from radiometer noise. It scales 
as B -1 / 2 for small bandwidths but steepens as larger bandwidths include frequencies where the flux density is higher. The 
difference between the total error and the standard error for the two-parameter fit is caused by the scattering term, tc- This 
example shows that a 3-parameter fit can remove much of the scattering when large bandwidths are used, allowing pulsars 
with large DMs to be included in a pulsar timing array Right: Similar plots for errors in DM when it is included in two and 
three-parameter fits. 
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Fig. 8. — Plot of standard and systematic errors for multifrequency fits to arrival times in the same format and for the 
same pulsar parameters as in Figure [7] except for a dispersion measure DM = 100 pc cm -3 and a maximum frequency of 
2 GHz. The left-hand panel is for a 100-m class telescope with >S S y S = 15 Jy and the right-hand panel is for an Arecibo-class 
telescope with Ssys = 4 Jy. 
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Fig. 9. — Plot of standard and systematic errors for multifrequency fits to arrival times in the same format and for the 
same pulsar parameters as in Figure \7\ except for a dispersion measure DM = 500 pc cm -3 and an upper frequency of 5 GHz. 
The left-hand panel is for a an Arecibo-class telescope with S B ys = 4 Jy) and the right-hand panel is for the full SKA with 
Ssys = 0.3 Jy. 



